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ABSTRACT 

Computer  simulation  techniques  were  used  to  determine 
equilibrium  positions  and  binding  energies  of  inert  gas  atoms 
implanted  in  a  tungsten  crystal  and  to  investigate  the  po- 
tential wells  around  these  equilibrium  positions  in  both  per- 
fect lattices  and  relaxed  lattices.   Stable  positions  were 
found  for  inert  gas  interstitials  near  lattice  atoms  in  the 
third  and  fourth  layers  of  the  crystal.   Interstitial  posi- 
tions near  atoms  in  the  first  and  second  layers  of  the  crys- 
tal appeared  to  be  unstable  if  they  exist  at  all.   As  a 
result  of  potential  well  studies,  it  was  concluded  that  the 
mechanism  associated  with  equilibrium  position  formation  was 
a  combination  of  local  liquefaction  of  the  lattice  structure 
and  interaction  of  the  interstitial  with  lattice  atoms. 
Equilibrium  positions  were  found  to  be  ill-defined  regions 
in  the  general  (110)  direction.   The  binding  energy  deter- 
mined for  an  interstitial  site  near  a  lattice  atom  in  the 
third  layer  of  the  crystal  was  in  excellent  agreement  with 
experimental  results. 
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I.   INTRODUCTION 

In  1968  Kornelsen  and  Sinha  [1]  of  the  National  Research 
Council  of  Canada  published  results  of  radiation-damage  ex- 
periments performed  on  tungsten.   In  these  experiments  a 
tungsten  surface  was  bombarded  with  ions  of  neon,  argon, 
krypton,  and  xenon  respectively.   Then,  while  the  tungsten 
was  heated,  gas  desorption  rates  were  measured  as  the  gas 
evolved  from  the  crystal.   The  resultant  desorption  spectrum 
was  interpreted  to  yield  a  binding  energy  spectrum  for  the 
trapped  particles.   In  1970  Professor  Don  E.  Harrison,  Jr., 
of  the  Naval  Postgraduate  School  undertook  the  modeling  of 
these  experiments  utilizing  computer  simulation  techniques 
in  order  to  provide  a  means  for  interpretation  of  their  re- 
sults.  Two  successive  thesis  research  efforts  [2,3]  have 
been  specifically  directed  toward  the  investigation  of  inert 
gas  implantation  in  a  tungsten  crystal.   It  was  anticipated 
that  a  corollary  to  the  successful  computer  simulation  of 
this  problem  might  possibly  be  an  improved  understanding  of 
the  interstitial  atom  stabilization  mechanism  in  tungsten, 
with  more  general  application  to  other  materials. 

The  investigations  reported  in  this  paper  are  a  con- 
tinuation of  work  begun  by  Vine  [21  and  Tankovich  [3]  under 
Harrison's  supervision.   The  simulation  procedures  followed 
were  a  combination  of  static  and  dynamic  approaches  to  the 
problem.   The  static  portion  of  the  problem  entailed 


interstitial  implantation  of  an  inert  gas  atom  in  a  tung- 
sten crystal  and  the  subsequent  relaxation  of  the  crystal 
until  the  equilibrium  position  of  the  interstitial  could  be 
ascertained.   In  the  dynamic  portion  of  the  problem  de- 
creasing amounts  of  energy  were  imparted  to  the  interstitial 
in  its  equilibrium  position  until  the  minimum  energy,  and 
direction,  which  still  allowed  the  interstitial  to  escape 
from  the  crystal  was  determined  (i.e.,  the  binding  energy 
of  the  particle) .   The  binding  energies  thus  determined 
provided  a  basis  for  comparison  with  the  results  of  Kornelsen 
and  Sinha. 
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II.   THE  NATURE  OF  THE  PROBLEM 

A.   THE  LATTICE  DYNAMICS  PROBLEM  IN  THE  COMPUTER 

For  a  little  over  a  decade  high  speed  digital  computers 
have  exhibited  their  usefulness  as  tools  for  investigation 
of  physical  systems.   Specifically,  the  inherent  periodicity 
and  order  of  crystal  lattices  have  made  the  study  of  lattice 
dynamics  particularly  adaptable  to  investigation  by  computer 
simulation.   It  is  a  relatively  simple  matter  to  "construct" 
a  crystal  of  the  desired  body-centered  cubic  or  face-centered 
cubic  structure  in  the  computer.   Various  types  of  point  de- 
fects can  also  be  "created"  in  the  crystal  with  relative 
ease.   A  vacancy  is  obtained  by  simply  removing  an  atom  from 
the  crystal,  while  interstitials  can  be  created  by  implanting 
an  additional  atom  in  the  crystal.   Two  types  of  interstitials 
have  been  used  in  investigations  of  lattice  dynamics,  atoms 
identical  to  the  crystal  lattice  atoms  ("self-interstitials" ) 
and  atoms  different  from  the  crystal  lattice  atoms 
(interstitials) . 

Although  a  crystal  lattice  containing  a  point  defect  can 
be  easily  represented  in  a  computer,  modeling  of  the  dynamics, 
which  allows  alterations  ^of  the  crystal  structure  resulting 
from  the  presence  of  the  point  defect,  is  more  involved. 
The  dynamic  portion  of  the  problem  of  computer  modeling  of 
lattices  can  be  characterized  by  four  key  decisions  which 
must  be  made. 


•First,  it  is  necessary  to  find  some  mathematical  re- 
lationship to  govern  the  interaction  among  the  atoms  of  the 
crystal.   This  interaction  is  a  complex  many  body  problem 
which  is  nearly  always  approximated  in  computer  simulations 
as  a  sum  of  appropriate  two  body  interactions.   With  this 
approximation  it  is  then  possible  to  represent  these  two 
body  interactions  by  some  type  of  potential  function,  which, 
in  turn,  can  be  used  to  determine  forces  on  individual  atoms. 

Secondly,  since  the  number  of  calculations  required  is 
directly  related  to  the  number  of  atoms  in  the  crystal,  a 
judicious  choice  of  crystal  size  must  be  made.   Large  crystals 
would  give  more  accurate  results,  but  would  also  require  more 
computational  time. 

A  third  key  problem  which  must  be  solved  results  from 
the  inability  of  a  digital  computer  to  perform  a  direct  in- 
tegration.  Since  the  lattice  dynamics  problem  is  most  often 
directed  at  a  determination  of  lattice  atom  positions  after 
some  type  of  interaction,  an  integration  of  the  equations  of 
motion  is  required.   A  choice  of  the  numerical  integration 
technique  to  be  used  must  therefore  be  made. 

Last,  and  probably  most  important,  since  an  iterative 
process  is  used  to  integrate  the  equations  of  motion,  the 
length  of  the  time  interval  of  the  iteration  must  be  chosen 
so  that  the  force  variations  within  this  time  interval  are 
small  and,  consequently,  stability  of  the  system  is  insured. 
At  first  glance,  this  suggests  that  only  very  small  timesteps 
should  be  taken,  but  this  would  require  excessive  computer 
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time.   Consequently,  some  variational  method  of  timestep 
determination  would  be  optimum. 

B.   HISTORICAL  DEVELOPMENT  OF  COMPUTER  SIMULATION  OF 
LATTICE  DYNAMICS 

The  four  basic  problems  of  computer  simulation  of  lattice 
dynamics  have  been  solved  in  various  ways  in  previous 
investigations . 

In  1960  Gibson,  Goland,  Milgram,  and  Vineyard  [4] 
(referred  to  hereafter  as  Gibson,  et.al)  of  the  Brookhaven 
National  Laboratory  published  the  first  complete  statement 
of  computer  simulation  procedure  as  a  tool  for  investigation 
of  radiation  damage  of  crystal  materials.   Due  to  the  com- 
plexity of  the  radiation  damage  problem  and  the  inadequacy 
of  analytical  methods  as  a  means  of  analysis  of  the  damage 
processes,  Gibson,  et.al.  turned  to  a  numerical  integration 
of  the  problem  utilizing  high  speed  computers.   Their  pio- 
neering work  gives  insight  into  such  pertinent  subjects  as 
crystal  size,  choice  of  potential  functions,  computational 
methods  for  solving  the  equations  of  motion,  and  timestep 
determination . 

Specifically,  their  research  utilized  copper  as  the 
crystal  material  since  its  relatively  simple  face-centered 
cubic  structure  and  its  widespread  use  in  actual  experimen- 
tation made  it  particularly  appropriate  for  initial  investi- 
gation.  A  Born-Mayer  repulsive  potential  was  used  to  describe 
atom-atom  interaction,  and  a  constant  cohesive  force  was  ap- 
plied to  each  atom  on  the  crystal  surface  to  balance  the 
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Born-Mayer  repulsion.   A  central  difference  iterative 
procedure  was  used  to  integrate  the  equations  of  motion. 
The  procedure  used  for  timestep  determination  is  of  parti- 
cular significance  since  the  principles  involved  are  still 
the  governing  criteria  for  choice  of  timestep  duration.   The 
fundamental  observation  was  made  that  the  greatest  stress 
on  the  crystal  system  was  a  result  of  the  strongest  atom- 
atom  interaction  in  the  system.   This  interaction  was  there- 
fore chosen  as  the  basis  for  timestep  determination.   Simply, 
the  timestep  duration  was  chosen  to  be  inversly  proportional 
to  the  velocity  created  by  the  strongest  atom-atom 
interaction. 

In  addition  to  verifying  the  applicability  of  computer 
simulation  techniques  to  radiation  damage  studies  and  pro- 
viding specific  information  on  collision  chains  and  focusing 
phenomena  in  crystallites,  the  work  of  Gibson,  et.al.[4] 
determined  that  the  only  stable  configuration  for  self- 
interstitials  in  a  face-centered  cubic  crystal  was  the  (lOO) 
split  configuration.   The  split  configuration  implies  that 
the  interstitial  causes  a  lattice  atom  to  split  away  from 
its  normal  lattice  position,  and  then  both  atoms,  the  inter- 
stitial and  the  lattice  atom,  share  (or  split)  the  normal 
lattice  site  when  equilibrium  is  reached. 

Johnson  and  Brown  [5]  confirmed  the  split  configuration 
as  the  only  stable  interstitial  position  in  their  studies 
of  copper  utilizing  a  Born-Mayer  repulsive  force  between 
nearest  neighbors  and  an  elastic  continuum  containing  the 
remainder  of  the  atoms  of  the  crystal.   Erginsoy,  Vineyard, 
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and  Englert  [6]  (referred  to  hereafter  as  Erginsoy,  et.al.) 
extended  these  calculations  to  the  body-centered  cubic  case 
by  performing  investigations  using  «  iron.   Computational 
techniques  paralleled  those  of  the  Brookhaven  group  except 
for  the  choice  of  a  potential  function.   After  experimen- 
tation with  a  Born -Mayer  potential  and  a  Morse  potential 
with  parameters  derived  by  Girifalco  and  Weizer  [7] , 
Erginsoy,  et.al.  [6]  settled  on  a  combination  potential  which 
utilized  an  exponentially  screened  Coulomb  potential  for 
close  approaches,  a  Born-Mayer  potential  for  first  nearest 
neighbor,    and  a  Morse  or  modified  Morse  potential  for 
higher  order  neighbors.   This  research  verified  the  split 
configuration  as  the  only  stable  configuration  for  body- 
centered  cubic  structures,  but  the  orientation  of  the  split 
was  found  to  be  in  a  (110)  plane,  (vice  a  (100)  plane  as  in 
the  face-centered  cubic  case) .   R.A.  Johnson  [8]  confirmed 
the  split  interstitial  orientation  in  similar  work  with  cc 
iron,  vanadium  and  tungsten. 

In  their  work  on  collisions  between  a  copper  atom  and 
a  copper  lattice,  Gay  and  Harrison  [9]  introduced  the  average 
force  procedure  for  integrating  the  equations  of  motion.   A 
complete  description  of  this  procedure  has  been  reported  by 
Gay,  Effron,  and  Harrison  [10] .   In  two  more  recent  efforts 
Johnson  and  Wilson  [11,12]  determined  new  potential  functions 
for  use  in  face-centered  cubic  and  body-centered  cubic  defect 
calculations  and  published  results  of  defect  calculations  of 
helium  in  various  metals. 
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Parallel  to  the  computer  simulation  efforts  described 
above,  Kornelsen  [13,14]  and  Kornelsen  and  Sinha  [1,15] 
have  published  the  results  of  extensive  experimentation 
involving  the  interaction  of  inert  gas  ions  with  tungsten. 
Under  Harrison's  supervision,  Vine  [2]  initiated  attempts 
to  devise  a  computer  model  of  Kornelsen' s  experiments  on 
neon  implantation  in  tungsten.   He  attempted  to  develop  a 
relationship  between  equilibrium  positions  of  tungsten 
interstitials  in  a  tungsten  crystal  determined,  first,  by 
using  a  Born-Mayer  repulsive  potential  to  describe  tungsten 
interstitial  interaction  with  a  tungsten  lattice  then,  by 
using  the  same  tungsten-tungsten  composite  potential  which 
was  assumed  between  atomsof  the  tungsten  crystal.   This  re- 
lationship would  then  have  been  applied  to  neon  equilibrium 
positions  derived  from  a  repulsive  potential  to  provide  a 
comparison  with  Kornelsen* s  data.   These  efforts  met  with 
little  success. 

Follow-on  work  by  Tankovich  [3]  utilized  the  static/ 
dynamic  approach  described  in  more  detail  in  Section  III-B. 
The  potential  function  used  in  these  investigations  was  a 
composite  Born-Mayer  and  Morse  potential  joined  by  the  best 
cubic  fit  in  the  area  of  intersection  which  had  previously 
been  developed  by  Harrison  and  Moore  [16] .   Static  program 
runs  confirmed  the  (110)  split  interstitial  for  helium, 
neon,  argon,  krypton,  and  xenon  point  defects  in  a  tungsten 
crystal  for  possible  interstitial   sites  in  planes  three 
through  six  of  the  ten  plane  crystal  used.   Preliminary 
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dynamic  testing  of  an  argon  equilibrium  position  in  plane 
four  of  the  crystal  was  begun  with  limited  success. 

Of  particular  significance  in  Tankovich's  work  [3]  was 
the  introduction  of  a  timestep  decrementing  process  into  the 
static  program.   In  previous  investigations  the  timestep 
duration  was  chosen  at  the  outset  of  the  problem  and  remained 
constant  through  all  computations.   At  this  author's  sug- 
gestion, a  timestep  decrementing  process  was  devised  which 
allowed  a  more  rapid  approach  to  the  equilibrium  positions 
with  a  concomitant  saving  of  computer  time  required  for 
computation.   (See  Section  III-A-3  for  a  discussion  of  the 
timestep. ) 
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III.   THE  SIMULATION  PROCEDURE 

A.   THE  MODEL 

1.   The  Crystal 

Two  factors  are  of  primary  importance  in  determining 
the  size  of  the  crystal  to  be  used  in  computer  simulation 
investigations.   First,  the  crystal  must  be  large  enough  to 
provide  realistic  results,  and,  second,  the  crystal  must  be 
as  small  as  possible  in  order  to  minimize  computer  time  re- 
quired for  calculation  of  atom-atom  interactions. 

After  experimenting  with  crystals  of  various  sizes, 
Tankovich  [3]  determined  that  a  crystal  with  ten  planes  of 
atoms  in  each  coordinate  direction,  which  is  equivalent  to 
five  unit  cells,  (10  x  10  x  10)  was  suitable  for  static  in- 
vestigations with  tungsten.   The  same  crystal  was  consequently 
used  in  the  investigations  reported  in  this  paper.   The  lat- 
tice unit,  or  distance  between  adjacent  (100)  planes  of 
atoms,  for  a  body-centered  cubic  tungsten  crystal  is  1.58A. 
All  distances  in  the  computer  simulation  were  measured  in 
lattice  units.   The  lattice  constant  for  the  tungsten  crystal 
is  3.16A  (or  two  lattice  units),  and  the  nearest  neighbor 
distance  is  ^  3   lattice  units. 

The  same  numbering  sequence  for  atoms  of  the  crystal 
that  was  employed  by  Tankovich  [3]  was  used  in  these  investi- 
gation.  (See  Figure  1.)   Atom  number  one  was  always  assigned 
to  the  interstitial  atom.   A  rectangular  coordinate  system 
was  placed  on  the  upper,  left  hand,  front  face  of  the  crystal. 
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The  positive  "x"  direction  was  chosen  to  the  right  from  the 
origin,  the  positive  "y"  direction  was  down,  and  the  positive 
"z"  direction  was  to  the  rear.  Atoms  were  numbered  consecu- 
tively, beginning  with  atom  number  2  at  the  origin  and  con- 
tinuing until  all  atoms  of  the  Y  =  0  plane  had  been  numbered. 
This  same  numbering  sequence  was  then  followed  for  each  Y 
plane  of  the  crystal  until  all  250  atoms  of  the  crystal  had 
been  numbered.   Figure  1  shows  the  numbering  of  the  atoms 
in  the  Y  =  0  and  Y  =  1  planes.   The  numbering  of  atoms  in  the 
remainder  of  the  planes  follows  the  same  pattern. 

For  the  dynamic  program  it  was  felt  that  a  smaller 
crystal  could  be  used  and  still  yield  meaningful  results. 
The  rationale  for  this  determination  was  based  on  the  go-no 
go  (i.e.,  escape  or  no  escape)  character  of  the  dynamic  pro- 
gram.  The  dynamic  program  essentially  provides  an  answer  to 
this  question  -  Will  an  interstitial  atom  escape  from  the 
crystal  if  it  is  given  a  specific  energy  and  directed  in  a 
specific  direction?   If  the  energy  dissipation  mechanism 
provided  by  collisions  with,  and  close  approaches  to,  the 
lattice  atoms  along  the  path  traveled  is  great  enough,  the 
interstitial  will  remain  in  the  crystal.   If  these  energy 
dissipation  mechanisms  are  not  large  enough  to  overcome  the 
kinetic  energy  of  the  interstitial,  the  interstitial  will 
escape.   Since  the  atoms  of  the  crystal  along  this  line  of 
motion  must  provide  the  mechanism  to  dissipate  the  inter- 
stitial 's  kinetic  energy,  for  escape  to  be  prevented,  only 
atoms  in  the  horizontal  planes  above  the  interstitial  and 
atoms  in  vertical  planes  within  a  few  lattice  units  of  the 
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interstitial  will  affect  the  possible  escape.   The  remainder 
of  the  atoms  of  the  crystal  would  not  have  time  to  react 
with  the  interstitial  or  affect  its  movement.   As  a  minimum, 
the  smaller  crystal  could  be  used  to  eliminate  excessively 
high  or  low  energies  from  consideration  at  a  considerable 
savings  of  computer  time  and,  instead,  give  a  limited  range 
of  energies  to  be  checked  by  dynamic  runs  using  the  entire 
crystal. 

To  implement  this  procedure  SUBROUTINE  PLUCK  was 
developed.   (See  Appendix  A  for  a  complete  discussion  of 
SUBROUTINE  PLUCK.)   Basically,  SUBROUTINE  PLUCK  uses  the  re- 
sults of  the  static  program,  but  causes  a  crystal  to  be 
printed  out  that  only  contains  the  interstitial  and  all  lat- 
tice atoms  from  two  planes  below  the  interstitial  to  the 
surface  plane  of  the  crystal  and  all  atoms  in  vertical  planes 
which  are  within  two  (or  three)  lattice  units  of  the  vertical 
plane  containing  the  lattice  site  shared  by  the  split  inter- 
stitial.  (See  Figure  2.)   The  savings  in  computer  time  re- 
sultingfrom  the  use  of  this  smaller  crystal  is  demonstrated 
by  the  following  comparison.   A  thirty  timestep  dynamic  run 
with  the  entire  10  x  10  x  10  crystal  (250  atoms)  used  approxi- 
mately three  and  one  half  minutes  of  computer  time.   A  thirty 
timestep  dynamic  run  using  a  7  x  5  x  7  crystal  (60  atoms) 
used  slightly  less  than  one  minute  of  computer  time  -  a  71% 
savings  in  computer  time! 

Most  of  the  dynamic  runs  made" during  the  course  of 
these  investigations  utilized  the  7x5x7  "PLUCK"  crystal. 

Final  confirmation  of  minimum  energies  was  determined  using 

the  entire  crystal. 
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2.   The  Potential  Function 

As  mentioned  previously,  the  many-body  interaction 
which  characterizes  actual  lattice  dynamics  is  approximated 
in  computers  by  many  two-body  interactions.   These  two-body 
interactions  are  represented  in  the  computer  by  some  type  of 
central,  pairwise  potential  function.   Various  types  and 
combinations  of  potential  functions  have  been  considered  for 
use  in  computer  simulation  investigations  of  lattice  dynamics. 

The  choice  of  the  potential  function  must  be  made 
with  due  consideration  to  the  range  of  applicability  of  the 
potential  function,  the  correlation  of  potential  function 
parameters  with  observable  properties  of  the  material  being 
investigated,  and  the  amount  of  computer  time  required  for 
calculations  using  that  potential  function.   The  potential 
function  used  in  these  investigations  was  the  composite 
Born-Mayer  potential  and  Girifalco  and  Weizer  Morse  potential 
used  by  Tankovich  [3,16]  .   This  composite  potential  is  con- 
structed as  follows: 

a.  Region  1- (r  <  1 .  5A) 

The  atom-atom  interaction  at  close  approach  is 

represented  by  a  Born  Mayer  repulsive  potential  of  the  form, 

$  .   =  exp  (A+B  r .  .  )  (1) 

13      ^v      13 

where  $. .  is  the  interaction  energy  between  particles  i  and 

j  and  r. .  is  the  distance  between  particles  i  and  j. 

b.  Region  2-(1.5A"  <  r  <  2.0&) 

This  portion  of  the  potential  function  is  ob- 
tained by  computing  the  best  cubic  equation  between  the  value 
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of  the  Born  Mayer  potential  at  1.5&  and  the  value  of  the 

Morse  Potential  at  2.0A. 

c.   Region  3-(2.0A"  <  r  <  5.38&) 

For  equilibrium  and  greater  separations,  a 

Morse  potential  of  the  form, 

|±j  =  D[exp{-2a(rij-ro)}-  2  exp{-a(rij-ro) }] 


(2) 


where  $. .  is  the  interaction  energy  between  particles  i  and 
j,  D  is  the  dissociation  energy  of  the  particles  i  and  j,  r._ 
is  the  distance  between  particles  i  and  j,  and  r   is  the 
equilibrium  separation,  is  used.   The  Morse  potential  was 
computed  so  that  the  tail  of  the  function  was  truncated  to 
zero  at  5.38A.   This  effectively  meant  that  atoms  out  to  the 
fourth  nearest  neighbor  were  included  in  calculations  of 
interaction  energies.   Girifalco  and  Weizer  [7]  had  pre- 
viously determined  Morse  potential  parameters  for  tungsten 
and  other  elements  which  included  interactions  out  to  the 
150th  nearest  neighbor.   Use  of  the  complete  function, 
however,  would  have  required  an  excessive  amount  of  computer 
time  for  calculation.   Additionally,  contributions  to  the 
interaction  energy  of  all  atoms  beyond  the  fourth  nearest 
neighbor  is  essentially  insignificant  for  our  calculations. 
Many  computer  simulations  of  lattice  dynamics  of  body- 
centered  cubic  materials  have  utilized  potential  functions 
which  only  included  first  and  second  nearest  neighbors  with 
satisfactory  results  [6,11].   To  check  the  adequacy  of  this 
potential  function  in  describing  this  crystal  system,  the 
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largest  binding  energy  observed  in  the  crystal  model 
(-8.283  eV)  may  be  compared  with  the  experimentally  de- 
termined heat  of  sublimation  for  tungsten  (-8.8  eV)   reported 
by  Harrison  and  Magnuson  [17].   This  agreement  within  5.9 
percent  was  considered  satisfactory  for  these  simulations. 
3.   The  Timestep 

a.   General  Discussion  of  the  Timestep 

Most  of  the  early  simulations  of  lattice  dynamics 
used  a  central  difference  method  as  the  numerical  procedure 
of  integrating  the  equations  of  motion  of  the  atoms  in  the 
crystal.   See  Gibson,  et.al  [4]  or  Gay,  Effron  and  Harrison[10] 
for  an  explanation  of  the  central  difference  method  of  numer- 
ical integration.   The  investigations  reported  in  this  paper, 
however,  used  the  average  force  method  which  is  completely 
described  by  Gay,  Effron,  and  Harrison  [10] .   Inherent  to 
both  numerical  methods  of  integration  is  the  replacement  of 
time  derivatives  in  the  equations  of  motion  with  a  finite 
time  difference,  i.e.,  the  timestep  interval.   As  mentioned 
previously,  the  strongest  interatomic  interaction  places  the 
greatest  demand  on  the  system;  consequently,  the  timestep 
duration  is  usually  determined  through  consideration  of  the 
strongest  interatomic  interaction  of  the  system.   These  in- 
vestigations followed  the  procedure  of  Gay,  Effron,  and 
Harrison  [10]  and  determined  the  timestep  duration  by  a  con- 
sideration of  the  maximum  displacement  that  the  most  ener- 
getic atom  of  the  crystal  was  allowed  to  move.   This  parameter, 
referred  to  hereafter  and  in  all  computer  programs  as  DTI,  is 
determined  as  follows: 
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(1)  The  equation  of  motion  of  i  atom  of  the 
crystal  during  timestep  interval  AT  can  be  written  in  the 
form 

Xi(t+AT)=  xi(t)  +  [vi(t)+  <Fi>AT/2m]AT  (3) 

or  in  the  equivalent  form 

Axi  =  (vi+  <Fi>AT/2m) AT.  (4) 

(2)  Rearranging  terms  of  equation  (4)  yields 

AT  =  Ax±/(vi+  <F.)AT/2m)  (5) 

where  AT  is  the  timestep  interval,  Ax.  is  the  displacement 

4-Vi 

of  the  i    atom  during  the  timestep,  v.  is  the  velocity  of 
the  i    atom,  and  (F . )  is  the  average  force  on  the  i    atom. 

From  equation  (5)  it  can  be  seen  that  the 
timestep  interval  is  a  function  of  both  the  kinetic  energy 
and  the  force.   Rather  than  solving  this  quadratic  equation 
for  the  timestep  interval,  the  average  force  method  considers 
separately  the  cases  when  energies  dominate  forces 
(v.  »<F .  )AT/2m)  and  when  forces  dominate  energies 
(<F.)AT/2m  »  v.).   Solutions  for  each  of  these  cases  yields 
a  different  value  for  the  (next)  timestep  interval. 

In  energy  dominant  cases, 

AT  =  Ax. /v.  (6) 

or,  expressed  in  the  form  of  energies, 

AT  =  Axi(m/2Tm)"2  (6a) 

and,  in  force  dominant  cases, 

AT  =     (2m    Axi/<Fi»3s,  (7) 
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where  Tm  of  equation  (6A)  represents  the  kinetic  energy 
of  the  atom  of  the  crystal  with  the  greatest  kinetic  energy, 
and  (F^)  of  equation  (7)  represents  the  average  force  on 
the  atom  with  the  maximum  force,  the  Ax.  of  each  equation 
becomes  DTI. 

In  these  investigations,  DTI  was  set  at  the 
beginning  of  the  program.   Ideally,  the  proper  choice  of 
timestep  duration  for  the  first  timestep  and  the  proper 
choice  of  DTI  would  lead  to  a  smooth  movement  of  the  inter- 
stitial to  its  equilibrium  position  in  a  manner  similar  to 
the  movement  of  a  critically  damped  oscillator. 

It  was  anticipated  that  energies  would 
dominate  forces  in  early  timesteps  which  would  lead  to 
timestep  determination  from  the  energy  equation,  equation 
(6A) .   At  some  point  in  the  crystal  relaxation  procedure, 
energies  should  have  been  dissipated  to  the  point  where 
further  timestep  determination  would  become  force  dependent 
(equation  (7) ) . 

b.   The  Average  Force  Method  and  Timestep  Determination 
In  the  average  force  method  of  integration  of  the 
equations  of  motion,  velocities  (i.e.,  energies)  of,  and 
forces  among,  all  atoms  of  the  crystal  are  computed  with  the 
atoms  in  their  initial  positions.   Based  on  these  forces  and 
energies  and  the  initial  timestep  duration,  new  positions 
for  all  atoms  of  the  crystal  are  computed.   The  forces  at  the 
new  positions  are  then  averaged  with  the  forces  at  the  ori- 
ginal positions  to  determine  the  average  force,  and  hence  the 
final  positions  of  all  atoms  at  the  end  of  the  first  timestep. 
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In  the  meantime,  the  maximum  force  exerted  on  an 
atom  of  the  crystal  in  its  original  position  and  the  maximum 
force  exerted  on  an  atom  of  the  crystal  in  its  final  (aver- 
aged) position  are  used  with  the  DTI  and  equation  (7)  to 
calculate  two  possible  alternatives  for  the  next  timestep 
interval.   Likewise,  the  energies  of  the  most  energetic  atoms 
in  both  original  and  final  positions  are  used  with  DTI  and 
equation  (6A)  to  calculate  two  other  possible  alternatives 
for  the  next  timestep  interval.   These  four  alternatives  are 
then  compared,  and  the  smallest  is  chosen  as  the  next  time- 
step  interval. 

c.   Procedures  Used  to  Determine  DTI 

Vine[2]  utilized  a  constant  DTI  in  all  of  his 
computations.   Since  the  movement  of  an  interstitial  to  an 
equilibrium  position  cannot,  in  general,  be  characterized  by 
a  small  range  of  energies  and  forces  and  since  DTI  should  be 
closely  correlated  with  energies  and  forces  at  least  in  ap- 
propriate regions,  the  use  of  a  constant  DTI  in  all  calcu- 
lations made  the  initial  choice  of  DTI  extremely  critical. 
Success  could  only  be  attained  by  resorting  to  small  DTI ' s 
and  concomitant  excessive  program  run  times  (>  100  timesteps)  , 

Tankovich  [3]  obtained  more  satisfactory  results 
by  successively  decrementing  DTI  for  each  timestep.   This 
procedure  allowed  long  timesteps  in  early  portions  of  a  run 
when  the  interstitial  was  far  from  its  equilibrium  position. 
As  equilibrium  was  approached,  the  decrementing  process  had 
progressed  to  the  point  such  that  significantly  smaller  and 
smaller  timesteps  were  taken  allowing  a  smooth  arrival  at  the 
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equilibrium  position.   Additionally,  this  was  accomplished 
at  a  considerable  savings  of  computer  time  (~  30  timesteps) . 

Although  this  decrementing  process  for  DTI  pro- 
vided considerable  improvement,  on  occasion,  the  final  few 
timesteps  used  such  a  small  DTI  that  practically  no  atom 
movement  was  discernible.   During  these  investigations  this 
situation  was  alleviated  by  incorporating  a  minimum  DTI  into 
the  decrementing  process.   This  insured  that  atom  movement 
was  still  discernible  near  equilibrium  and  assisted  in 
guarding  against  possible  false  assumption  of  equilibrium 
because  of  the  relatively  small  movement  observed  under  the 
continuous  decrementing  process. 

In  an  attempt  to  more  fully  understand  the 
mechanisms  of  the  static  solution,  the  computer  program  was 
adjusted  so  that  the  maximum  force,  the  atom  upon  which  this 
force  was  exerted,  and  the  four  "new"  possible  timesteps  for 
each  timestep  interval  were  printed  out  after  each  run.   It 
was  observed  that  the  timestep  calculations  based  upon  "new" 
and  "old"  energies  were  overwhelmingly  the  basis  for  timestep 
determination.   To  insure  that  the  force  dependence  of  the 
timestep  was  not  being  unduly  disregarded,  the  program  was 
adjusted  so  that  DTI  was  determined  solely  as  a  function  of' 
the  minimum  of  the  two  forces.   These  "force  calculations" 
of  equilibrium  positions  agreed  with  "energy  calculations" 
within  0.03  lattice  unit. 
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B.   THE  PROGRAMS 

Although  the  basic  computational  procedures  contained 
in  the  computer  programs  for  the  static  and  dynamic  portions 
of  the  problem  were  essentially  the  same,  an  understanding 
of  the  subtle  differences  between  the  programs  and  an  ap- 
preciation for  several  computational  tools  and  procedures 
used  in  the  programs  have  to  be  gained  prior  to  further  dis- 
cussion of  the  actual  investigation  and  results. 

1.   The  Static  Program 

In  the  static  program,  a  tungsten  crystal  of 
appropriate  size  was  created  in  the  computer.   An  inter- 
stitial atom  was  then  implanted  at  a  chosen  site  within 
the  crystal.   Potential  energies  and  mutual  forces  of  all 
atoms  were  computed.   The  crystal  was  then  allowed  to  relax 
in  appropriately  chosen  timesteps.   At  the  end  of  each  time- 
step  an  energy  dissipation  mechanism  was  introduced  in  the 
form  of  a  predetermined  damping  factor  which  was  used  to 
decrease  all  velocity  components  of  the  atoms  of  the  crystal. 
The  next  timestep  interval  was  then  computed,  and  the  process 
was  repeated  until  an  indicated  number  of  timesteps  had  been 
completed.   If  an  equilibrium  position  had  been  reached, 
positions  and  energies  of  all  atoms  in  the  crystal,  including 
the  interstitial,  and  other  pertinent  data  could  be  punched 
out  on  cards  for  later  use  in  the  dynamic  program. 
a.   Equilibrium  Positions  of  Interstitials 

In  interpreting  the  results  of  the  static  simu- 
lations it  was  necessary  to  arrive  at  some  criterion  to  use 
as  a  determining  factor  for  the  interstitial ' s  final 
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equilibrium  position.   It  was  expected  that  the  crystal 
would  relax  around  the  split  interstitial  site  forming  a 
"pocket"  within  the  crystal  which  would  interrupt  the 
periodicity  of  the  lattice.   Stability  in  this  configuration 
was  determined  to  have  been  reached  when  the  atoms  of  the 
deformed  crystal,  including  the  interstitial,  had  been 
allowed  to  relax  (i.e.,  adjust  to  the  presence  of  the 
interstitial)  to  the  point  where  their  kinetic  energies 
were  all  below   thermal.   (<  0.025  eV) 

If  more  than  one  possible  equilibrium  position 
met  this  criterion,  it  was  felt  that  different  positions 
within  the  "potential  well"  of  the  equilibrium  site  were 
probably  being  observed.   The  position  in  which  the  inter- 
stitial atom  had  the  smallest  amount  of  potential  energy  was 
then  chosen  as  the  equilibrium  position  for  the  lattice  site 
under  investigation. 

b.   Handling  of  Oscillations  Near  Equilibrium 

While  performing  the  "force  calculations",  some 
runs  began  with  the  interstitial  moving  toward  an  apparent 
equilibrium  until  an  oscillation,  or  "rattle",  developed 
about  the  suspected  equilibrium  position  in  the  "z"  direction 
It  was  first  confirmed  that  this  "rattle"  was  solely  in  the 
"z"  direction  (i.e.,  no  significant  "x"  or  "y"  displacement) 
by  extending  the  computation  from  30  to  200  timesteps.   It 
was  shown  that  in  a  (100)  direction  in  a  body-centered  cubic 
a  narrow  potential  energy  minimum  was  observed  at  (100) 
planes.   The  region  of  the  (110 )  line  between  the  octahedral 
void  and  the  reference  lattice  site  was  contained  within 
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this  potential  energy  minimum.   With  this  information  as 
justification,  it  was  determined  that  the  velocity  of  the 
interstitial  in  the  "z"  direction  could  be  completely  damped 
at  the  end  of  each  timestep  to  hasten  the  determination  of 
the  true  equilibrium  position.   This  technique  was  useful 
in  restricting  computer  run  time  whenever  an  obvious  "rattle" 
in  the  "z"  direction  developed. 
2.   The  Dynamic  Program 

The  initial  step  of  the  dynamic  program  was  the  re- 
creation in  the  computer  of  the  tungsten  crystal,  including 
the  interstitial,  after  relaxation  of  the  crystal  had  taken 
place  and  the  interstitial  had  come  to  an  equilibrium  posi- 
tion. This  was  accomplished  by  utilizing  the  output  of  the 
static  program  as  the  input  to  the  dynamic  program. 

Before  continuing  with  the  problem,  it  was  necessary 
to  realize  that  the  minimum  energy  required  for  the  inter- 
stitial to  escape  from  the  crystal  would  be  a  function  of 
the  path  of  escape.   To  account  for  this  "direction  depend- 
ence", impact  points  were  chosen  in  the  surface  plane  of  the 
crystal  perpendicular  to  the  direction  of  escape.   The  nega- 
tive "y"  direction  was  always  used  as  the  direction  of 
escape.   Energy  was  then  imparted  to  the  interstitial  which 
was,  in  turn,  "aimed"  at  a  specific  impact  point.   The  inter- 
stitial was  then  allowed  to  travel  through  the  crystal  using 
a  timestep  procedure  based  on  a  constant  DTI.   (The  timestep 
procedure  in  the  dynamic  program  is  simpler  than  in  the 
static  program  since  energies  are  dominant  throughout.)   Each 
impact  point  (i.e.,  direction  of  escape)  was  subsequently 
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tested  in  the  same  manner  using  the  same  initial  energy. 
If  the  interstitial  escaped  from  the  crystal  in  any  di- 
rection tested,  it  was  assumed  that  the  initial  energy 
applied  to  the  interstitial  was  greater  than  the  binding 
energy  for  that  particular  ion  in  that  particular  equili- 
brium position.   In  this  case,  the  initial  energy  was  de- 
creased and  another  survey  of  the  impact  points  was  taken. 
When  the  minimum  initial  energy  which  still  allowed  the 
interstitial  to  escape  was  determined,  the  binding  energy 
for  that  particular  ion  in  that  particular  equilibrium 
position  was  known. 

3.   Equilibrium  Sites  Viewed  As  Potential  Wells 

One  means  of  modeling  equilibrium  positions  of 
interstitials  in  a  crystal  lattice  is  to  consider  each  pos- 
sible equilibrium  position  as  a  three  dimensional  potential 
well.   A  foreign  atom  migrating  through  the  lattice  which 
reaches  the  confines  of  one  of  these  potential  wells  with 
sufficiently  small  energy  would  fall  into  the  potential  well 
and  remain  there.   In  conjunction  with  this  research,  the 
character  of  these  potential  wells  was  also  investigated. 

Two  different  approaches  were  taken  to  investigate 
the  potential  well  aspect  of  the  problem.   First,  the  static 
program  was  used  to  investigate  behavior  of  interstitials 
implanted  at  various  positions  around  a  previously  deter- 
mined equilibrium  position  in  a  perfect  lattice  to  determine 
which  of  these  positions  sought  equilibrium.   The  second 
method  offset  the  interstitial  from  its  equilibrium  position 
in  the  relaxed  crystal  and  then  allowed  the  crystal  to  undergo 
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further  relaxation  to  determine  whether  the  offset  inter- 
stitial would  seek  the  original  equilibrium  position. 
4.   Determination  Of  Possible  Interstitial  Sites 

In  a  body-centered  cubic  material  there  are  twelve 
possible  locations  for  a  (110)  split  interstitial.   In  an 
infinite  crystal  these  sites  are  equivalent;  that  is,  no  one 
site  can  be  distinguished  from  another.   When  a  finite  crystal 
is  considered,  the  presence  of  a  crystal  surface  allows 
identification  of  three  distinct  interstitial  sites.   An  "A" 
site  is  located  in  a  (110)  plane  which  contains  the  lattice 
site  about  which  the  split  occurs  and  is  closer  to  the 
crystal  surface  than  the  shared  lattice  site.   A  "C"  site 
is  also  in  a  (110)  plane  containing  the  shared  lattice  site, 
but  is  located  below  the  shared  lattice  site.   A  "B"  site 
is  located  in  the  same  (100)  plane  parallel  to  the  surface 
that  contains  the  shared  lattice  site. 

a.   Implantation  Procedure 

In  the  static  program  the  interstitial  was 
initially  implanted  in  an  offset  position  in  the  direction 
of  the  expected  equilibrium  position.   This  procedure  al- 
lowed a  smooth  movement  toward  the  suspected  equilibrium 
position  in  a  minimum  number  of  timesteps.   For  the  heavier 
interstitial  atoms  (xenon  and  tungsten)  the  lattice  atom 
was  also  offset  in  its  direction  of  suspected  movement  as 
a  result  of  the  presence  of  the  large  interstitial. 
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IV.   PRESENTATION  OF  DATA 

A.   STATIC  SIMULATION 

1.  Initial  Simulations 

After  the  minimum  DTI  procedure  had  been  incorporated 
into  the  DTI  decrementing  process,  investigations  of  tungsten 
interstitials  in  a  tungsten  crystal  were  made.   In  particular, 
a  comparison  between  tungsten  movement  under  the  influence 
of  an  attractive  potential  and  tungsten  movement  under  the 
influence  of  a  repulsive  potential  was  sought.   Equilibrium 
positions  determined  separately  with  these  two  potentials 
gave  agreement  within  0.02  lattice  unit. 

Investigations  were  then  extended  to  verify  equili- 
brium positions  near  lattice  sites  89  and  64  which  had  pre- 
viously been  obtained  by  Tankovich  [3].   Positions  for  argon, 
neon,  krypton,  and  xenon  were  examined  for  each  site,  and 
results  agreed  with  Tankovich' s  data  within  0.02  lattice 
unit,  see  Table  I.   As  the  mass  of  the  interstitial  atom 
decreased,  its  equilibrium  position  moved  from  the  shared 
lattice  site  location  previously  reported  [7,8]  toward  the 
center  of  the  octahedral  void. 

2.  Site  39 

Investigations  were  then  directed  toward  determining 
equilibrium  positions  of  a  split  interstitial  near  lattice 
site  39,  which  is  one  layer  below  the  surface  of  the  crystal. 
Calculations  were  made  separately  for  argon,  neon,  krypton, 
and  xenon  interstitials.   At  the  end  of  thirty  timesteps, 
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TABLE  I 
Interstitial  and  Lattice  Atom  Displacements  from  Reference  Site 


Interstitial 

Interstitial 
Site 

Interstitial 
Displacement 

Lattice  Atom 
Displacement 

Ax 

Ay 

Ax 

Ay 

Neon 
Argon 
Krypton 
Xenon 

89C 
89C 
89C 
89C 

-0.86 
-0.80 
-0.77 
-0.48 

+0.87 
+0.80 
+  0.77 
+0.48 

+0.06 
+0.19 
+0.24 
+0.41 

-0.06 
-0.19 
-0.24 
-0.42 

Tungsten 

89C 

-0.38 

+0.39 

+0.35 

-0.35 
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all  atoms  of  the  crystal  in  each  case  had  kinetic  energies 
below  thermal,  which  indicated  that  an  equilibrium  had  been 
reached.   An  observation  was  made,  however,  that  tended  to 
abrogate  this  determination.   The  lighter  interstitials, 
neon  and  argon,  exhibited  a  definite  affinity  for  the  sur- 
face of  the  crystal.   In  all  three  possible  sites  (A,B,  and 
C) ,  argon  and  neon  interstitials  were  found  to  seek  posi- 
tions significantly  (~  0.1  lattice  unit)  closer  to  the  sur- 
face than  the  expected  (110)  split.   Final  positions  for 
argon  and  neon  interstitials  in  site  39A  were  located  less 
than  three  tenths  of  a  lattice  unit  below  the  surface  of 
the  crystal.   Similar  behavior,  but  to  a  lesser  degree, 
was  also  observed  with  the  heavier  interstitials,  krypton 
and  xenon.   Additionally,  although  the  kinetic  energy  cri- 
teria indicated  that  an  equilibrium  had  been  reached,  the 
small  velocities  that  were  available  at  the  end  of  each 
timestep  were  in  such  a  direction  as  to  allow  escape  from 
the  crystal  if  these  velocities  could  be  maintained.   In 
fact,  although  velocities  were  halved  at  the  end  of  each 
timestep,  a  comparison  of  velocities  over  the  last  five 
timesteps  of  the  computer  run  showed  that  the  negative  "y" 
(direction  of  escape)  velocity  component  of  the  interstitial 
at  the  end  of  timestep  thirty  was  actually  only  twenty  per- 
cent lower  than  the  same  velocity  component  at  the  end  of 
timestep  twenty  five.   In  other  words,  even  with  fifty  per- 
cent damping  applied  during  each  timestep,  velocities  actu- 
ally decreased  by  only  twenty  percent  over  _5_  timesteps.   It 
was  also  observed  that  during  these  last  five  timesteps  the 
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maximum  movement  allowed  by  the  most  active  atom  of  the 
crystal  (the  DTI)  had  reached  its  minimum  value,  0.0005 
lattice  unit.   These  observations  suggested  that  the  small 
movement  allowed  during  these  timesteps  combined  with  the 
damping  factor  might  be  "forcing"  the  interstitial  to  ex- 
hibit equilibrium  criteria.   It  was  concluded  that  equili- 
brium positions  near  site  39  were  unstable,  if  they  exist 
at  all. 

3.   Site  14C 

Until  these  investigations  little  thought  was  given 
to  the  possibility  of  an  equilibrium  position  in  site  14C. 
With  no  precise  knowledge  available  concerning  the  actual 
quantitative  value  of  the  damping  experienced  by  a  foreign 
interstitial  implanted  in  a  crystal  lattice  and  with  the 
possibility  of  equilibrium  sites  near  lattice  site  39, 
some  credence  had  to  be  given  to  the  possibility  of  equili- 
brium positions  near  surface  atoms  of  the  crystal.   It  was, 
therefore,  decided  to  investigate  the  possibility  of  an 
equilibrium  position  in  site  14C.   It  was  postulated  that 
in  the  case  of  light  atoms  (neon,  for  example)  the  inter- 
stitial position  for  the  14C  site  should  be  deeper  in  the 
crystal  than  the  39A  site,  as  a  result  of  the  greater  re- 
pulsion of  the  lattice  atom  with  which  the  site  was  shared. 
Simulations  showed  that  the  14C  interstitial  site  was  located 
at  a  distance  of  0.43  lattice  unit  below  the  surface  of  the 
crystal  while  the  39A  interstitial  site  was  located  at  a 
distance  of  0.26  lattice  unit  below  the  crystal  surface. 
When  the  displacement  of  this  interstitial  site  was  compared 
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with  the  64C  interstitial  site  (0.84  lattice  unit  below 
the  next  higher  lattice  plane) ,  the  effect  of  the  surface 
of  the  crystal  on  interstitials  in  close  proximity  to  the 
surface  was  clearly  exhibited.   Velocity  characteristics 
similar  to  those  of  site   39A  were  also  observed  indicating 
that  the  amount  of  damping  applied  and  the  DTI  could  be 
forcing  the  interstitial  to  exhibit  equilibrium  properties 
in  this  site.   In  general,  it  can  be  said  that  interstitial 
equilibrium  sites  in  the  first  two  layers  of  the  crystal 
are  ill  defined  in  position  if  they  exist  at  all. 
4.   Force  Calculations 

As  mentioned  previously,  in  the  course  of  these 
investigations  it  was  determined  that  energies  nearly  al- 
ways dominated  forces  for  the  timestep  ranges  used  in  the 
calculations,  and,  consequently,  the  timestep  was  nearly 
always  chosen  as  a  function  of  the  energy  of  the  most  ener- 
getic atom.   To  ascertain  whether  the  use  of  timesteps  deter- 
mined by  maximum  forces  would  yield  better  or  significantly 
different  results,  the  static  program  was  modified  so  that 
the  timesteps  were  determined  strictly  as  a  function  of  the 
maximum  force.   This  was  accomplished  by  first  printing  out 
the  maximum  force  observed  during  each  timestep.   By  sur- 
veying these  maximum  forces,  several  values  of  force  were 
chosen  to  be   used  as  test  forces.   As  the  program  was  sub- 
sequently run,  the  maximum  force  in  each  timestep  was  compared 
with  these  test  forces,  and  a  DTI  for  that  timestep  was  as- 
signed based  on  the  results  of  that  comparison.   This 
assigned  value  of  DTI  was  then  used  to  compute  the  next 
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timestep  interval.   These  calculations  of  equilibrium 
positions  based  on  forces  agreed  with  previous  energy  cal- 
culations of  equilibrium  positions  within  0.03  lattice 
unit. 

It  should  be  realized  here  that  the  terminology 
"force  calculations"  and  "energy  calculations"  do  not  imply 
any  significant  difference  in  method  for  determining  equili- 
brium positions.   They  are  merely  two  different  procedures 
for  determining  the  maximum  displacement  which  will  be  al- 
lowed by  -the  most  energetic  atom  of  the  crystal  during  the 
next  timestep.   The  importance  of  this  parameter  (DTI)  is 
in  the  effect  it  has  in  insuring  the  smooth  movement  of  the 
interstitial  to  its  equilibrium  position. 
5.   Investigation  of  Potential  Wells 

Viewing  interstitial  equilibrium  positions  as  po- 
tential wells  of  varying  depths  is  a  convenient  means  of 
modeling  the  entrapment  of  foreign  atoms  by  lattice  struc- 
tures.  In  order  to  minimize  the  effects  of  the  crystal 
surface  on  investigations,  interstitial  site  89C,  which  had 
exhibited  good  stability  and  almost  perfect  (110)  splitting 
in  previous  testing,  was  chosen  for  investigation  of  the 
characteristics  of  interstitial  potential  wells. 

a.   Potential  Well  Studies  in  a  Perfect  Lattice 
The  first  approach  to  the  study  of  the  inter- 
stitial potential  well  utilized  essentially  the  same  method 
that  had  been  originally  employed  to  locate  equilibrium  posi- 
tions.  It  was  postulated  that  any  interstitial  atom  that  was 
implanted  in  a  perfect  lattice  at  or  near  the  coordinates  of 
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the  equilibrium  position  which  had  been  previously  deter- 
mined for  interstitial  site  89C  would  seek  the  same  equili- 
brium site.   Interstitial  atoms  could  be  implanted  further 
and  further  from  the  previously  determined  equilibrium 
position  until  they  no  longer  returned  to  it.   In  this  man- 
ner a  map  of  the  size  of  the  potential  well  around  site  89C 
could  be  obtained. 

The  results  obtained  in  investigations  of  xenon 
in  site  89C  are  presented  in  Figure  3.  In  this  figure  each 
coordinate  intersection  represents  an  implantation  site 
which  was  tested.  The  arrows  drawn  from  these  implantation 
sites  indicate  the  direction  of  movement  of  the  interstitial 
from  that  specific  implantation  site.  The  tip  of  the  arrow 
represents  the  interstitial  position  after  thirty  timesteps. 

In  considering  the  data  obtained  in  these  in- 
vestigations, two  significant  observations  were  made.   First, 
an  interstitial  implanted  at  the  previously  determined  equili- 
brium position  (coordinates  (4.52,  3.48)  in  Figure  3)    moved 
from  this  position  to  another  position  which  also  met  the 
criteria  for  equilibrium.   Secondly,  all  other  implantation 
sites  also  moved  to  positions  which  met  equilibrium  criteria. 
An  analysis  of  the  program  print  out  data  provided  informa- 
tion of  secondary  importance.   All  implanted  atoms  exhibited 
an  initial  movement  in  the  general  (110 >  direction  even 
though  final  positions  were  not  necessarily  in  that  direction, 
Similar  investigations  were  conducted  using  argon  and  neon 
interstitials  with  similar  results. 
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In  light  of  this  unexpected  behavior  in  the 
perfect  lattice  the  stability  of  a  foreign  atom,  inserted 
into  the  crystal  as  a  replacement  for  the  lattice  atom  in 
site  89,  was  examined.   Investigations  were  made  using  neon, 
argon,  krypton,  and  xenon  as  the  replacement  atoms.   The 
data  from  these  runs  are  tabulated  below. 


TABLE  II 


Results  of  Replacement  Runs  in  Site  89 


Replacement 
Atom 

Neon 
Argon 
Krypton 
Xenon 


Y  Displacement  (L.U.) 
After  30  Timesteps 

-0.049 
-0.036 
-0.029 
-0.020 


Final  Kinetic 
Energy  (eV) 

0.0000 
0.0008 
0.0001 
0.0000 


Final  Potential 
Energy  (eV) 

1.0444 
3.5118 
3.5120 
1.3297 


The  y  displacement  values  are  sufficiently  small 
that  site  89  must  be  presumed  to  be  a  stable  replacement  site, 
This  conclusion  is  further  confirmed  by  observing  that  the 
potential  energies  are  also  significantly  lower  than  those 
obtained  for  the  same  atoms  when  acting  as  interstitials, 
see  Table  III. 


TABLE  III 
Interstitial  Potential  Energies  for  Site  89C 


Interstitial 
Atom 

Neon 
Argon 
Krypton 
Xenon 


Final  Potential 
Energy  (eV) 

4.7 
16.3 
15.5 

9.1 
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To  further  check  the  stability  of  this  site, 
a  kryptron  replacement  atom  was  then  initially  offset  from 
the  site  89  coordinates,  0.7  lattice  unit  in  the  "x"  and 
"y"  directions,  and  the  program  was  run  again.   The  krypton 
atom  moved  precisely  along  the  (110)  line  back  toward  site 
89.   At  the  end  of  thiry  timesteps,  the  krypton  atom  was 
located  0.4  lattice  unit  in  the  "x"  and  "y"  directions, 
"x"  and  "y"  velocity  components  were  still  in  the  direction 
of  site  89,  and  potential  energy  had  decreased  from  54  eV 
in  the  initial  offset  position  to  8.7  eV  at  the  end  of 
thirty  timesteps. 

(1)  Interpretation  of  Results 

Interpretation  of  the  results  reported  in 
the  previous  section  led  to  a  re-examination  of  the  concept 
of  the  equilibrium  positions  of  interstitials  and  raised  a 
question  concerning  the  validity  of  solely  using  relaxed 
crystals  for  determinations  of  equilibria. 

The  equilibrium  positions  obtained  by  using 
a  perfect  crystal  were  a  function  of  the  implantation  site. 
This  dependency  of  the  final  equilibrium  position  on  the 
implantation  site  suggested  that  the  movement  of  the  inter- 
stitial to  an  equilibrium  site  can  not  be  explained  in  terms 
of  the  lattice  structure  "forcing"  the  interstitial  to  its 
lowest  energy  configuration.   Rather,  these  results  sug- 
gested that  the  actual  mechanism  is  a  process  of  local 
"liquefaction"  of  the  lattice  as  it  adjusts  to  the  presence 
of  the  interstitial  combined  with  interstitial  movement  caused 
by  interaction  with  the  lattice  atoms.   This,  in  turn, 
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suggested  that  a  more  accurate  means  of  localizing  equili- 
brium positions  might  be  to  utilize  a  crystal  which  had 
already  been  allowed  to  "relax",  or  adjust  to  the  presence 
of  the  interstitial.   (See  next  section.) 

The  movement  of  atoms  from  nine  different 
implantation  sites  in  an  area  two  tenths  of  a  lattice  unit 
per  side  to  nine  different  positions  which  exhibit  equili- 
brium criteria  indicated  that  equilibrium  positions  are  not 
precise  positions  coordinate  -  wise  within  the  time  range 
of  these  computations. 

b.   Potential  Well  Studies  in  a  Relaxed  Lattice 

The  conclusions  reached  as  a  result  of  studies 
of  potential  wells  in  a  perfect  crystal  prompted  similar 
studies  in  a  crystal  which  had  been  allowed  to  adjust  to 
the  presence  of  the  interstitial.   An  argon  interstitial 
and  interstitial  site  89C  were  chosen  for  investigation. 
The  relaxed  lattice  was  obtained  by  using  the  results  of  the 
static  simulation  of  argon  in  site  89C  which  had  first  in- 
dicated that  an  equilibrium  position  had  been  reached.   These 
results  were  read  into  the  computer  as  initial  positions  for 
the  lattice  atoms  and  the  interstitial.   The  interstitial 
was  then  offset  from  its  position  and  a  new  static  simu- 
lation was  performed.   An  array  of  sites  was  chosen  for  in- 
vestigation in  this  manner.   The  results  of  these  simulations 
with  various  interstitial  offsets  are  shown  in  Figure  4. 
The  representation  in  Figure  4  is  analogous  to  the  repre- 
sentation in  Figure  3.   The  start  point  for  each  run  is 
numbered  (1-9  and  A,B). 
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The  interstitial  movement  depicted  in  Figure  4 
was  encouraging  in  many  respects.   An  argon  interstitial 
placed  successively  at  each  initial  site  (numbered  1-9  in 
Figure  4)  appeared  to  head  for  roughly  the  same  region  of 
the  crystal.   The  potential  energy  of  the  interstitial  at 
the  conclusion  of  each  run  ranged  in  value  from  8.8  to  9.8 
eV.   A  comparison  of  these  energies  with  the  energy  of  the 
interstitial  in  the  initial  equilibrium  position  found  in 
the  perfect  lattice  (~16  eV)  indicated  that  this  equilibrium 
region  is  much  more  stable  than  the  position  determined  pre- 
viously.  While  not  precisely  defined  in  position,  this 
equilibrium  region  was  located  in  roughly  the  (110)  di- 
rection from  the  shared  lattice  site. 

Since  this  equilibrium  region  appeared  to  be 
located  nearer  to  the  shared  lattice  site  (site  89)  than 
offset  start  points  1-9,  two  additional  computer  runs  (labeled 
A  and  B  in  Figure  4)  were  made  with  start  points  bordering 
this  equilibrium  region.   These  runs  indicated  that,  even 
in  a  relaxed  crystal,  the  equilibrium  reached  was  still 
somewhat  a  function  of  the  initial  position  of  the  inter- 
stitial.  In  both  of  these  runs,  equilibrium  was  reached  by 
movement  in  the  (110)  direction,  and  the  potential  energy  of 
the  interstitial  after  thirty  timesteps  (~  10  eV)  was  lower 
than  the  potential  energy  originally  determined  in  the  per- 
fect lattice. 

(1)  Interpretation  of  Results 

The  results  obtained  in  the  relaxed  lattice 
seemed  to  reinforce  the  conclusions  drawn  during  the  study 
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of  potential  wells  in  the  perfect  lattice.   Basically,  the 
relaxed  crystal  studies  gave  a  picture  of  equilibrium  posi- 
tions much  more  in  tune  with  what  might  reasonably  be  ex- 
pected in  nature.   It  seems  clear  that  the  equilibrium  site 
of  an  interstitial  is  an  ill-defined  region  in  the  general 
(110)  direction  from  the  shared  lattice  site.   Determination 
of  this  equilibrium  region  cannot  be  made  based  solely  on 
the  kinetic  energy  of  the  interstitial;  potential  energies 
must  also  be  considered. 

This  equilibrium  region  could  have  b een 
postulated  from  a  consideration  of  the  rms  displacement  of 
a  tungsten  atom  in  a  tungsten  lattice.   Houska  [18]  has 
measured  this  rms  displacement  to  be  0.049  lattice  unit  at 
300  K.   The  ordinary  thermal  activity  of  the  lattice  atoms, 
then,  can  be  expected  to  cause  the  equilibrium  positions  of 
interstitials  to  fluctuate  over  some  region.   This  equili- 
brium region  might  best  be  described  as  roughly  a  cylinder 
whose  axis  is  in  the  (110)  direction  and  whose  height  and 
radius  are  some  function  of  the  relationship  between  inter- 
stitial mass,  lattice  atom  mass,  and  the  interatomic  poten- 
tial function.   The  equilibrium  positions  determined  in  these 
investigations  appeared  to  be  centered  in  the  (110)  direction 
at  the  intersection  of  the  (110)  and  (100)  planes  through 
the  lattice  atom  and  covered  a  section  of  the  (110)  line 
on  the  order  of  0.2  of  a  lattice  unit  long.   This  seems  to 
be  a  reasonable  range  of  equilibrium  positions  for  the  re- 
latively light  argon  interstitial  in  a  tungsten  crystal. 
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The  results  of  these  investigations  have 
indicated  the  usefulness  of  conducting  implantation  studies 
in  relaxed  crystals.  The  movement  of  interstitials  to  equi- 
librium in  the  relaxed  crystal  was  much  less  dependent  upon 
initial  positions,  and  all  sites  investigated  showed  a  pre- 
ference for  positions  in  the  (110)  direction.  Such  was  not 
the  case  in  perfect  lattice  studies. 

B.   DYNAMIC  SIMULATIONS 

Concurrent  with  the  investigations  into  the  potential 
wells  of  the  equilibrium  sites,  dynamic  simulations  were 
conducted  using  the  argon  interstitial  in  site  64A.   Since 
Tankovich  [3]  had  investigated  an  argon  interstitial  in 
site  89C  and  had  determined  that  the  binding  energy  of   this 
site  was  in  the  range  2-4  eV,  4  eV  was  chosen  as  the  initial 
energy  for  the  dynamic  tests.   It  was  determined  that  the 
binding  energy  of  the  argon  interstitial  in  site  64A  was 
between  0.02  and  0.04  eV.   At  0.04  eV  escape  of  the  inter- 
stitial was  noted  for  several  of  the  twelve  impact  points 
tested.   At  0.02  eV  the  interstitial  first  moved  slightly 
toward  the  surface  of  the  crystal,  then  changed  directions 
and  moved  to  a  position  deep  inside  the  crystal  while  gaining 
considerable  energy.   This  phenomena  had  been  seen  and  in- 
terpreted previously  as  the  expected  behavior  of  a  stable 
site  in  which  the  interstitial  has  been  required  to  move 
too  large  a  distance  in  one  timestep.   Since  the  timestep  is 
constant  (0.1  lattice  unit)  in  the  dynamic  simulation,  an 
interstitial  that  is  oscillating  in  a  stable  configuration 
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can  be  required  to  move  out  of  its  stable  position.   When 
this  movement  is  directed  into  the  crystal,  it  can  be  assumed 
that  the  original  position  was  stable  . 

This  binding  energy  falls  in  the  range  of  the  first 
desorption  peak  measured  by  Kornelsen  and  Sinha  (see  Figure 
5  in  Ref .  1) .   This  result  indicated  that  interstitial  sites 
formed  with  lattice  atoms  in  the  third  crystal  plane  (cor- 
responding to  the  y  =  2  plane  of  these  calculations) ,  such 
that  the  interstitial  site  was  in  the  (110)  direction  above 
the  lattice  atom,  are  the  sites  nearest  the  surface  of  the 
crystal  stable  enough  to  entrap  measurable  amounts  of  argon. 
This  excellent  concurrence  with  experiment  further  sub- 
stantiated the  conclusion  drawn  earlier  in  these  investi- 
gations about  the  doubtful  character  of  positions  near  sites 
39  and  14  which  had  exhibited  equilibrium  criteria. 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

During  the  course  of  these  investigations  the  possibility 
that  the  DTI  decrementing  procedure  and  the  procedure  for 
incorporating  damping  into  the  problem  might  be  forcing  in- 
terstitials  to  exhibit  equilibrium  properties  arose.   At 
first,  the  DTI  was  thought  to  possibly  be  restricting  atom 
movement  to  such  an  extent  during  the  last  few  timesteps  of 
a  calculation  (i.e.,  when  DTI  had  reached  its  minimum  value 
-0.0005  lattice  unit)  that  atom  movement  and,  consequently, 
kinetic  energy  could  no  longer  be  used  as  equilibrium  cri- 
teria.  This  became  more  significant  when  the  small  DTI  and 
the  damping  factor  were  considered  together.   The  damping 
factor  seemed  particularly  appropriate  for  scrutiny  near  the 
surface  of  the  crystal  and  when  DTI  was  small,  since  the  pro- 
bability of  collisions  and  close  approaches  and  the  subsequent 
damping  of  atom  motion  could  be  expected  to  decrease  in  both 
instances.   The  excellent  agreement  with  experiment  of  the 
results  of  the  dynamic  and  static  simulations  suggests,  how- 
ever, that  the  DTI  and  damping  procedures  employed  might  be 
adequate  for  these  simulations  except  when  atom  movement  is 
particularly  close  to  the  surface  of  the  crystal.   It  is 
recommended  that  future  investigations  explore  the  use  of  a 
damping  factor  which  varies  with  decreasing  DTI  and  movement 
toward  the  surface  of  the  crystal. 
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These  investigations  have  shown  the  value  of  using  a 
relaxed  crystal  in  the  computer  simulation  of  lattice  dynamics. 
The  relaxed  crystal  simulation  indicated  that: 

(1)  The  equilibrium  position  initially  sought  by  an 
interstitial  is  a  function  of  implantation  position. 

(2)  The  mechanism  associated  with  the  establishment  of 
an  equilibrium  position  seems  to  be  a  combination  of  inter- 
atomic interaction  and  local  liquefaction  of  the  crystal 
structure  in  the  vicinity  of  the  interstitial. 

(3)  The  character  of  equilibrium  potential  wells  is  more 
readily  observable  in  a  relaxed  crystal. 

(4)  The  actual  equilibrium  position  of  an  interstitial 
seems  to  be  some  region  in  the  general  (110)  direction  from 
the  reference  (shared)  lattice  site. 

Again,  however,  for  determining  binding  energies,  the 
original  procedure  for  determining  equilibrium  positions  may 
yield  adequate  results.   Such  was  the  case  for  the  dynamic 
simulations  reported  here.   It  is  recommended  that  replace- 
ment interstitials  and  other  split  interstitial  sites  undergo 
dynamic  simulation  in  an  attempt  to  correlate  other  desorption 
peaks  as  determined  by  Kornelsen  and  Sinha  [1]  with  inter- 
stitial sites. 
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APPENDIX  A 
SUBROUTINE  PLUCK 

Subroutine  Pluck  was  developed  to  allow  the  use  of  a 
crystal  smaller  than  the  original  model  in  dynamic  testing. 
This  subroutine  was  developed  so  that  the  only  pieces  of 
information  required  as  input  to  the  subroutine  were  the 
number  of  crystal  planes  desired  in  the  "x"  and  "z" 
directions  and  the  number  of  the  lattice  site  desired  at 
the  center  of  the  PLUCK  crystal.   Once  the  size  of  the 
PLUCK  crystal  is  decided  upon,  the  subroutine  stores  the 
original  numbers  of  the  atoms  in  the  PLUCK  crystal  in  an 
array.   This  allows  reference  to  any  atom  by  its  original 
number  throughout  computation.   The  atoms  of  the  PLUCK 
crystal  were  numbered  consecutively,  and  the  number  of  atoms 
in  the  PLUCK  crystal  was  assigned  a  variable  name.   In  this 
manner,  a  minimum  of  adjustment  was  required  in  the  dynamic 
program  when  the  PLUCK  crystal  was  used. 

SUBROUTINE  PLUCK  is  included  in  this  appendix  in  its 
most  general  form.   Parameters  are  listed  below  for  two 
different  sizes  of  PLUCK  crystals.   Sizes  refer  to  the 
number  of  "x"  and  "z"  planes  in  the  PLUCK  crystal.   "Y"  planes 
from  the  surface  of  the  crystal  through  the  two  planes  below 
the  plane  of  the  lattice  site  under  investigation  are  always 
included. 
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PLUCK  PARAMETERS  FOR  DIFFERENT  SIZED  CRYSTALS 

(crystal  sizes  refer  to  number  of  planes  of  atoms  in  x  and 

z  directions) 

Parameter 

IXNEW 

IZNEW 

NM1 

Nil 

NX1 

NII31 

NII41 

NIINC1 

IIINC1 

NMINC1 

NXINC1 

NMINC2 

NX2 

NM2 

NI2 

NII32 

NII42 

NUNC  3 

IIINC2 

NMINC3 

NXINC2 

NMINC4 


7x7  Crystal 

5x5  Crystal 

7 

5 

7 

5 

5 

3 

8 

14 

16 

4 

10 

16 

5 

9 

1 

3 

4 

2 

4 

2 

9 

9 

-1 

1 

9 

9 

4 

4 

8 

8 

4 

10 

11 

15 

2 

2 

3 

3 

3 

3 

16 

4 

1 

-1 
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•SUBROUTINE    PLUCK' 


C0MM0N/C0M1/RX(500) , RY (5  00 ) , RZ ( 500 ) »LCUT(500), 
1LL,LD, I TYPF,NVAC 

C0MMQN/C0M4/IXt IY, IZ* SCX,SCY,SCZ, IDEEPtDlX,DlY,DlZt 
UVACXtl  VACY.IVACZ 

COMMON /COM  10/ I  XN FW , I YNEW , I ZNE W ♦ I  I 

COMMON/COM11/RXNEWI (250),RYNEWI ( 250 ) , RZNEW I ( 250) , 
1KEEP(250) ,NNUM(250) 
1500     IXNEW=IXNEW 

IYNEW=IVACY+3 

IYNEW=I YNEW 

NM=NM1 

NI=NU 

11  =  2 

MM=0 

NX=NX1 

NII3=NI 131 

NI  I4=NI 141 

IF( IYNEW.E0.3)     GO    TO    1514 

IF{  IYNEW.E0.5)     GO    TO    1514 
1505    DO    1509    1=1  I ,NM 

NNUM( I ) =NI 
1509    NI=NI+1 

NI=NI+NI INC1 

1 1 =11  +  1 1 1  NCI 

NM=NM+NMINC1 

IF     (II  .  LE.NX)     GO    TO    1505 

NX=NX+NXINC1 

NI=NI+NI 14 

MM=MM+1 

IF(IYNEW.EO.M«l)     GO    TO     1600 

NM=NM+NMINC2 

GO  TO  1515 

1514  NX=NX2 
NM=NM2 
NI=NI2 

NI  I3  =  NI 132 
NI  I4=NI 142 

1515  DO    1520    1=1  I ,NM 
NNUM(  I  )=NI 

1520    NI=NI+1 

NI=NI+NI INC3 

II=II+I I INC2 

NM=NM+NMINC3 

IF ( I  I .LE.NX)     30    TO    1515 

NX=NX+NXINC2 

NI=NI+NI 13 

MM=MM+1 

IF( IYNEW.EQ.MM )    GO    TO    1600 

NM=NM+NMINC4 

GO    TO    1505 
1600     11=11-1 

RXNEWI  (1)=RX(1  ) 

RYNEWI ( 1)=RY(1 ) 

RZNEWI  (1  )  =  RZ(1  ) 

KEEP(l) =1 

NNUM( 1)=1 
1700    DO    1750     1=2,11 

RXNEWI ( I )=RX( NNUM1 I ) ) 

RYNEWK  I  )  =  RY(NNUM(  I  )  ) 

RZNEWI (I  )  =  RZ(NNUM(  I  )  ) 

KEEP(  I  ) =NNUM(I  ) 
1750    CONTINUE 

RETURN 

END 
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APPENDIX  B 
COMPUTER  PROGRAM  GLOSSARY 
ALPHA:  .      Input  Morse  potential  parameter. 
BSAVE:       Target  Mass/  (target  mass  +  bullet  mass) ; 

distributes  potential  energy  between  target  and 

bullet. 
BIND:        Negative  of  the  total  potential  energy  (TPOT) 

at  time  zero. 
BMAS:        Mass  of  bullet  in  AMU. 

BULLET:      Alpha-numeric  array  for  point  defect  material. 
CFO,  CFl  ,  CF2:  Force  parameters  of  cubic  fit  between  Morse 

and  Born -Mayer  functions. 
CGB1,  CGB2:  Morse  potential  parameters. 
CGD1,  CGD2:  Morse  potential  parameters. 
CGF1,  CGF2:  Morse  force  parameters. 
CPO,  CPl,  CP2,  CP3:  Potential  parameters  of  cubic  fit 

between  Morse  and  Born-Mayer  functions. 

CVD:  CVR  X  10    :  Converts  lattice  units  to  meters. 

-19 
CVE:  1.6  X  10    :  Converts  electron  volts  to  Joules. 

CVED:  CVE/CVD:  A  ratio  used  to  avoid  repeated  division. 

-27 
CVM:  1.67  2  X  10    :  Converts  atomic  mass  units  to  kilograms, 

CVR:  LU  in  angstroms;  converts  lattice  units  to  angstrom 

units. 
DIX,  DIY,  DIZ:  Displacement  coordinates  for  location  of 

interstitial  from  reference  atom,  NVAC. 
DCON:        Input  Morse  potential  parameter. 
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DDTIF:       The  minimum  value  that  DTI  is  allowed  to  assume. 

DDTII:       The  initial  decrement  of  DTI. 

DFF:         ROE-DIST.  The  distance  closer  than  ROE  that  an 

atom  is  to  the  primary. 
DIST:        Distance  between  any  two  atoms. 
DLPE:  TLPE-TLPE0:  The  change  in  total  local  potential  energy 

since  time  zero. 
DRX,  DRY,  DRZ:  x,y,z  components  of  DIST. 
DT:  Length  of  a  timestep  in  seconds. 

DTE,  DTE1:   The  two  possible  alternatives  of  the  timestep 

computed  from  maximum  energies. 
DTF,  DTF1:   The  two  possible  alternatives  of  the  timestep 

computed  from  maximum  forces. 
DTI:         Number  of  lattice  units  most  energetic  atom  may 

move  in  one  timestep. 
DTI  (I),  DT2  (I),  DT3  (I),  DT4  (I):  Vector  arrays  which  save 

the  possible  choices  of  timestep  determined  in 

the  "energy"  method. 
DTSTEP  (I) :  Vector  array  which  saves  the  timestep  interval 

chosen  for  each  timestep. 
DTOD:  DT/CVD  —  A  ratio  used  to  avoid  repeated  division. 
DTOM:DT/PTMAS  —  A  ratio  used  to  avoid  repeated  division. 
DX(I),  DY(I),  DZ(I):  Change  in  position  of  ith  atom  from 

initial  position  at  time  zero. 
EMAX:        The  maximum  energy  encountered  in  any  cycle. 
EV:  Primary  energy  in  electron  volts. 

EVR:         Primary  energy  in  kiloelectron  volts. 
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EXAf  EXB:    Input  Born-Mayer  potential  function  parameters 

for  the  target. 
F2:  Square  of  the  force  on  a  specific  atom. 

FA:  The  component  force  increment  on  an  atom. 

FDTI :  DTI  X  CVD:  A  parameter  used  to  determine  DT  by 

maximum  energy  method. 
FM:  A  small  number  used  in  checking  potential  energy 

zero  point. 
FM2:  FM  squared. 
FMAX:        Maximum  total  force  on  the  most  stressed  atom 

in  the  crystal. 
FOD:  FORCE/DIST:  A  ratio  used  to  avoid  repeated  division. 
FORCE:       Numerical  value  of  the  force  function  with  a 

variable  parameter. 
FORMAX  (I):  Vector  array  which  saves  the  maximum  force  in 

each  timestep. 
FX(I),  FY(I)/  FZ(I):  x,y,z,  components  of  total  force  on 

an  atom. 
FAX:         Born-Mayer  force  function  parameter. 
HBMAS:  %   BMAS :  A  ratio  used  to  avoid  repeated  division. 
HDTOD:  \   DTOD :  A  ratio  used  to  avoid  repeated  division. 
HDTOM:  %  DTOM:  A  ratio  used  to  avoid  repeated  division. 
HDTOMB:  ^  DTOMB :  A  ratio  used  to  avoid  repeated  division. 
HTMAS:  \   TMAS :  A  ratio  used  to  avoid  repeated  division. 
II:         Variable  in  cubic  fit  subroutine. 
13:         Variable  in  cubic  fit  subroutine. 
IDEEP:       First  fixed  layer. 
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IFMXAM  (I) :  Vector  array  which  saves  the  atom  number  which 

experiences  the  maximum  force. 

Alpha  numeric  array  for  program  title. 

Alpha  numeric  array  for  Morse  function  parameters 

Alpha  numeric  array  for  bullet  element. 

Alpha  numeric  array  for  type  and  orientation  of 

crystal. 

Alpha  numeric  array  for  target  element. 

Number  of  atoms  in  a  crystal  using  subroutine 

PLUCK. 

Number  of  free  (mobile)  layers. 

Odd-even  integer  used  to  determine  atom  site 

establishment. 

Subscript  value  of  atom.   Void  in  subroutines 

STEP  and  ENERGY. 

Parameter  that  determines  whether  or  not  a  self 

defect  is  to  be  given  a  repulsive  potential  or 

a  composite  attractive  -  repulsive  potential. 

A  parameter  used  to  shut  down  the  program. 

Unsealed  fixed  point  x  coordinate  used  in 

lattice  generation. 

Odd-even  integer  used  to  determine  atom  site 

establishment. 

Parameter  used  to  determine  the  type  of  point 

defect:  vacancy,  interstitial,  or  replacement. 
IX,  IY,  IZ:  Number  of  x,y,z  planes  of  cyrstal. 
IXNEW,  IYNEW,  IZNEW:  Number  of  x,y,  and  z  planes  in  the 

PLUCK  crystal. 


IH1 
IH2 
IHB 
IHS 

IHT: 
II: 


I  LAY: 
IN: 

IP: 

IQ: 


ISHUT: 
IT: 

ITT: 

I TYPE: 
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J2:  Variable  in  the  cubic  fit  subroutine. 

JJ:  Parameter  in  the  BCC(lll)  lattice  generation 

subroutine. 
JT:         Unsealed  y  coordinate  used  in  crystal  generation. 
JTS:        Variable  used  to  establish  atom  sites. 
JTT:         Variable  used  to  establish  atom  sites. 
KEEP  (I) :    Vector  array  which  saves  the  original  atom 

numbers  of  the  PLUCK  crystal. 
KF:  Final  K  in  LOCAT  (K)  assigned  to  an  atom. 

KT:         Unsealed  z  coordinate  used  to  establish  atom 

site. 
LCUT  (I) :    Used  to  identify  an  ith  atom  which  is  not 

included  in  calculations. 
LD:  The  highest  numbered  atom  in  the  mobile  layers. 

LL:  The  highest  numbered  atom  in  the  entire  crystal. 

LOCAT  (K) :   Dimensioned  variable  that  remembers  the  numbers 

of  the  atoms  within  a  radius  ROEL  of  the  primary 

at  time  zero. 

Variable  associated  with  each  of  the  nine  lattice 

generator  subroutines. 

One  number  higher  than  the  order  of  the  fit 

between  the  Born -Mayer  and  Morse  potentials, 

always  4  in  this  simulation. 
ND:         Data  output  increment,  in  numbers  of  timesteps. 
NEW:         Parameter  used  to  determine  whether  or  not  atom 

numbers  have  been  stored  in  LOCAT  (K) . 
NNUM  (I) :    Vector  array  used  in  Subroutine  PLUCK  to  re- 
number atoms. 


LS: 


MCRO: 


54 


NPAGE 
NRUN: 

NS: 
NT: 
NTT: 
NVAC: 


PAC: 

PBMAS: 
PEXA,  PEXB 

PFPTC : 

PFXA: 

PKE  (I): 

PLANE: 

POT: 

PPE  (I): 

PPEINT  (I) 

PPENCK: 

PPENEG  (I) 

PPEPCK: 

PPEPOS  (I) 


Page  numbering  variable. 

Parameter  used  to  determine  whether  or  not  to 

read  additional  data  cards. 

Initial  print  statement  timestep  number. 

Timestep  number. 

Timestep  number  limit  before  shutdown. 

An  atom  number  used  to  establish  point  defects 

or  used  as  a  reference  point  for  interstitial 

placement. 

Parameter  for  bullet  force  function  correction. 

Primary  mass  in  kilograms. 

Input  Born-Mayer  potential  function  parameters 

for  the  bullet-target  interaction. 

Primary  force  function  evaluated  at  ROE. 

Primary  force  function  parameter. 

Kinetic  energy  of  the  ith  atom. 

Alpha-numeric  array  for  lattice  orientation. 

Potential  energy  between  two  atoms. 

Potential  energy  of  the  ith  atom. 

Vector  array  that  saves  the  difference  in 

potential  energy  before  and  after  implantation. 

Potential  energy  check  value  which  determines 

potential  energy  decreases  which  will  be  printed. 

Vector  array  which  saves  potential  energy 

decreases. 

Potential  energy  check  value  which  determines 

which  potential  energy  increases  will  be  printed, 

Vector  array  which  saves  potential  energy  in- 
creases. 
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PPESAV  (I) :  Vector  array  which  saves  the  initial  potential 
energy  of  lattice  atoms  before  implantation. 

PPKEEP  (I) :  Vector  array  which  saves  potential  energy  dif- 
ferences between  perfect  crystal  and  relaxed 
crystal. 

PPTC:        Primary  potential  function  evaluated  at  ROE. 

PTE  (I) :     Total  energy  of  the  ith  atom  (potential  + 
kinetic) . 

PTMAS:       Target  mass  in  kilograms. 

RE:  Input  Morse  potential  parameter. 

RO:  Spacing  constant  in  BCC (110)  lattice  generation 

subroutine. 

ROE:         Nearest  neighbor  distance. 

ROE  2:       ROE  squared. 

ROEA:        Maximum  cut-off  for  Born-Mayer  potential. 

ROEB:        Minimum  cut-off  for  Morse  potential. 

ROEC:'       Maximum  cut-off  for  Morse  potential. 

R0EC2:       ROEC  squared. 

ROEL:        Radius  inside  of  which  local  potential  energy 
is  found. 

R0EL2:       ROEL  squared. 

ROEM:        ROE-DTI,  Region  in  which  modification  of 
repulsive  force  must  be  made. 

RX(I)/  RY(I),  RZ(I):  x,y,z,  coordinates  of  an  ith  atom 
at  any  time. 

RXI(I),  RYK(I),  RXI(I):  x,y,z,  coordinates  of  an  ith  atom's 
initial  position. 
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RXK(I),  RYK(I),  RZK(I):  x,y,z  coordinates  of  temporary 

position  of  an  ith  atom  during  force  cycle. 
RXNEWI,  RYNEWI,  RZNEWI :  Vector  arrays  which  contain  the 

x,y,z,  coordinates  of  the  atoms  of  the  PLUCK 

crystal. 
RXSAVE,  RYSAVE,  RZSAVE:  x,y,z,  coordinates  of  the  impact 

point  in  the  dynamic  program . 
SAVE:        h   POT. 

SCX,  SCY,  SCZ:  x,y,z,  coordinate  scale  factors. 
SSCZ:    -    A  z  scale  factor  used  for  the  FCC  (HI)  lattice 

generator  subroutine. 
START:       An  optional  timing  variable,  not  used  in  this 

simulation. 
SUM:         Variable  in  cubic  fit  subroutine. 
TARGET:      Alpha-numeric  array  for  target  material. 
TSAVE:       Bullet  mass/(target  mass+bullet  mass); 

distributes  potential  energy  between  target 

and  bullet. 
TE:  Total  energy  of  all  crystal  atoms  (kinetic  + 

potential) . 
TEMP:        Temperature  of  lattice  in  degrees  K elvin  not 

used  in  this  simulation. 
TFAC:        A  time  factor  ratio  used  to  determine  DT  by 

maximum  force  method. 
TFACB:       TFAC  for  the  bullet. 

THERM:       Thermal  energy  of  atom  not  used  in  this  simu- 
lation. 
TIME:        Elapsed  problem  time  in  seconds. 
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TLPE:        Total  local  potential  energy  of  atoms  within 

a  radius  ROEL. 
TLPE0:       TLPE  at  time  zero. 
TMAS:        Target  atom  mass  in  AMU. 
TPKE:        Total  kinetic  energy  of  all  crystal  atoms. 
TPOT:        Total  potential  energy  of  all  crystal  atoms. 
VSS:         Storage  variable  for  velocity  components. 
VS  (1) ,  VY(1),  VZ(1):  x,y,z  components  of  ith  atoms  velocity. 
X,  Y,  Z:     Unsealed  coordinates  used  in  crystal  generation, 
XNVAC,  YNVAC,  ZNVAC :  The  initial  displacement  (in  LU)  of 

atom  NVAC. 
YLAX(I):     Relaxation  in  -y  direction  of  ith  layer  in  L.U. 
ZP:  Floating  point  form  of  JTT. 
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APPENDIX  C 
ERROR  IN  THE  LITERATURE 
It  was  discovered  during  the  course  of  these  investi- 
gations that  Equation  (20)  of  Ref.  10  (corresponding  to 
equation  6A  of  the  report)  was  incorrect. 
Equation  (19)  of  Ref.  10  is, 

AT  =  Axj/vi 

where  AT  is  the  timestep  interval,  ax.  is  the  displacement 
of  the  i   particle,  and  v.  is  the  velocity  of  the  i 
particle. 

To  express  the  timestep  in  terms  of  the  energy  of  the 
particle  with  the  maximum  energy,  Ax.  is  defined  as  the 
displacement  of  the  particle  with  the  greatest  energy  and 
is  replaced  by  the  symbol  D.   Tm  is  this  particle's  energy, 
which  is  the  largest  kinetic  energy  at  end  of  each  timestep. 

If  T   is  expressed  as 

m      r 

T   =  h   mv  , 
mm 

substitution  into  equation  (19)  yields, 

AT  =  D/(2  Tm/m)  2, 

°r     i  AT  =  D(m/2  Tm)"^  , 

vice  ^ 

AT  =  D(2m/T  )  2.   (Equation  (20) Ref.  10) 

It  should  be  noted  that  this  error  has  no  affect  on  the 

calculations  reported  here.   D  (DTI  in  these  calculations) 

is  a  parameter  which  has  been  specifically  chosen  in  order 
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to  provide  a  timestep  interval  over  which  there  are  no 
significant  force  changes.   The  factor  of  2  which  is  intro- 
duced in  this  correction  would  merely  have  required  a  sub- 
sequent alteration  of  the  DTI ' s  chosen  so  that  the  timestep 
interval  would  remain  essentially  the  same. 
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Figure  1.   Atom  Numbering  Procedure. 
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THE     FUNCTION    OF       SUBROUTINE       PLUCK 


PLUCK 
CRYSTAL 
(7x5x7) 


ORIGINAL 
CRYSTAL 
(10x10x10) 


Figure  2.   The  Function  of  Subroutine  Pluck 

62 


a: 
o 

UJ 

o 

> 

< 

< 

X 

_l 

LlI 
GQ 

< 

h~ 

Z 

z 

o 

u_ 

z 

or 

UJ 

UJ 

X 

a. 

-  o 

t  I 

CO  H 

O  UJ 

°-  3 


Figure  3.   Xenon  Behavior  in  a  Perfect  Lattice. 
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Figure  4.   Argon  Behavior  in  a  Relaxed  Lattice 
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COMPUTER  PROGRAM 
FOR  STATIC  SIMULATIONS 

STATIC  SIMULATIONS  WERE  USED  TO  INVESTIGATE  THE 
EOUILIBRIJM  POSITIONS  OF  INTERSTITIAL  ATOMS  IMPLANTED  IN 
A  TUNGSTEN  CRYSTAL.  FOUR  DIFFERENT  CONFIGURATIONS  OF  THE 
STATIC  PROGRAM  WERE  USED  AS  THE  INVESTIGATIONS  PROGRESSED. 
THESE  CONFIGURATIONS  WERE  USED  TO 

(1)  DETERMINE  EQUILIBRIUM  POSITIONS  BY  'ENERGY' 
CALCULATIONS, 

(2)  DETERMINE  EQUILIBRIUM  POSITIONS  BY  'FORCE' 
CALCULATIONS, 

(3)  PRINT  OUT  THE  SMALLAR  CRYSTAL  DETERMINED  BY 

SUBROUTINE  PLUCK  FOR  USE  IN  INITIAL  DYNAMIC 
TESTING,  AND 

(4)  INVESTIGATE  THE  POTENTIAL  WELLS  IN  THE  RELAXED 
CRYSTAL. 

THE  'ENERGY'  CONFIGURATION  OF  THE  STATIC  PROGRAM  IS  PRESENT- 
ED BELOW  WITH  DIFFERENCES  FROM  THIS  PROGRAM  REQUIRED  BY 
OTHER  CONFIGURATIONS  INCLUDED  AND  DISCUSSED  AT  APPROPRIATE 
POINTS.  ADDITIONALLY,  BRIEF  COMMENTS  ARE  INCLUDED  AT  VARIOUS 
POINTS  TO  CALL  ATTENTION  TO  SIGNIFICANT  PROCESSES  OF  THE 
PROGRAM. 

C 

C  DIMENSION    VECTOR    ARRAYS    USED    EXCLUSIVELY    IN    THE    MAIN 

C  PROGRAM.     THE    DIMENSIONING    SCHEME    FOR    EACH    CONFIG- 

C  URATION     IS    GIVEN    BELOW    IN    ITS    ENTIRETY. 

C 

C  'FORCE    CONFIGURATION' 

DIMENSION    VX(500),VY(500) , VZ1500) ,PKE(500) 

DIMENSION    DX(500) ,DY(500),DZ( 50  0),PTE(500) 

DIMENSION    PPESAV(500),PPEINT(500),PPEPOS(5  00), 
1PPENEG( 500) 

DIMENSION    FSTACC(IOO) ,FSSACC(100 ) , FORMA X { 100 ) , 
1IFMXAM( 100) 

DIMENSION    DTFSP(50 ) , DT I  ASP ( 50 ) , DTESP ( 50 ) , DTI VSP ( 50) , 
1DETK50J 

DIMENSION    PPKEEP(500) 

DIMENSION    DTSTEP(200) 


C  'OTHER    CONFIGURATIONS' 

DIMENSION    VX(500),VY(500),VZ(500),PKE(500) 

DIMENSION    DX(500) , DY(500),DZ(5OO),PTE(500) 

DIMENSION  PPESAV(500),PPEINT(500),PPEPOS(5  00), 
1PPENEG( 500) 

DIMENSION  FSTACCUOO) ,FSSACC(1D0 ) ,FORMAX( 100), 
1IFMXAM( 100) 

DIMENSION    DTI ( 50 ) , DT2 ( 50 ) , DT3 ( 50 ) , DT4{ 50) 

DIMENSION    PPKEEP(500) 

DIMENSION    DTSTEP(200) 
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c 

C  PRESCRIBE    COMMON    STORAGE    OF    VARIABLES    AND    VARIABLE 

C  ARRAYS    REQUIRED    IN    SUBROUTINES. 

C 

COMMON/ COM1/RX ( 500 ) ,RY ( 500 ) , RZ ( 500),LCUT{  5  00)  , 
1LL,LD, ITYPE,NVAC 

C0MM0N/C0M2/HK  20)  ,IH2C 8) , IHS( 10) , IHB(6)  ,IHT(6)  , 
1 TARGET (4 JtTMASt BULLET (4), BMAS , PL AN E, TEMP , THERM t 
1DDTI I»DDTIF 

C0MM0N/C0M3/RXI ( 500) , RYI ( 5  00) ,RZ I ( 500) ,C VR ,E VR , 
1NT,TIMF,DT,DTI  ,  I  LAY , RXK ( 50 0 ) , RYK ( 500 ) ,RZK( 500) 

C0MM0N/CCM4/IX, IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z, 
1IVACX. IVACY, IVACZ 

COMMON/CCM5/ROE,ROE2,ROEM,EXA,EXB,PEXA,  PEXB, FXA,PFXA, 
1IQ,TSAVE,3SAVE 

COMMON/ C0M6/FXI 5  00 ),FY(500),FZ(50  0),PAC,PFPTC,FM 

COMMON/COM7/PPTC,TP0T,PPE(1000) , T LP E , ROEL , R0EL2 , N EW 

COMM0N/COM8/ROEA,ROE3,R0EC,ROEC2,CP0,CPl ,CP2,CP3, 
1CF0,CF1 ,CF2,CGD1,CGD2,CGB1, CGB2 , CGF1 , CGF2 

CCMMCN/C0M9/XNVACYNVACZNVAC 

COMMON/COM10/IXNEW, IYNEWt IZNEW, II 

C0MM0N/C0M11/RXNEWI (25  0) ,RYNEWI (2  50 ) , RZNEW I ( 250 ) , 
1KEEP(250) tNNUM(250) 

COMMON/COMA/     A(4,5),MCR0 

C 

C  LIST    FORMAT    STATEMENTS    OF    ALL    READ    COMMANDS. 

C 

9010  FORMAT  (20A4) 

9020  FORMAT( 8A4,3F8. 5,2F5.2) 

90  30  F0RMAT(4A4,3F8.5,6A4,F6.2) 

9040  F0RMAT(F6.3t5X,  15 » 5  14  ,  4X , 3 F5 , 3  I  2 ) 

9041  F0RMAT(2|r6.4,6F6.3) 

9042  FORMAT (F8.4,F8.4 ) 

9  0  50    FORMAT(  1 0A4 , A 4 , 41 3 , F 8. 4 , 1  4 , F 5 . 4 ) 
9052     FORMAT     (6(F6.3)) 

C 

C      LIST  FORMAT  STATEMENTS  OF  ALL  WRITE  COMMANDS. 

C 

9610  FORMAT (1H1) 

9620    F0RMAT(47X« ' SUMMARY    OF    ATOMS' // ,35X , 8A4 , ' ,    NT    =«I4,//, 

13C     ATOM  POSITION  BIND    ENERGY  •),//) 

963  0    F0RMAT(3(I 5t3F6.2,F8.4»8X)) 
9640    FGRMAT(/4X,F10.3,25H    EV, TOTAL    KINETIC    ENERGY , ,F1 0 .3 , 

127H    EV, TOTAL    POTENTIAL     ENERGY, F  10.3,  •     EV,     REDUCTION1,/ 

1/60X, 'RADIUS    =' ,F5.2, i 
96  50     FORMAT! 10  5X,4HPAGE, 13, /, 1H1) 
9660    FORMAT (/       •  ATOM  DX  DY  DZ 

1VX  VY  VZ  KE  PE  TE',/) 

9670    F0RMAT(1I8,3F10.3,3F10.1,3F10.4    ) 
9680    FORMATC     SHARP     DT     DECR  E  AS  E  '  ,  2  E10  .3  ) 

9690  FORMATl 14, 3F5.2, 14) 

9691  F0RMAT{9F8.4) 

9692  FORMAT( IX, 14,/) 

9693  F0RMAT(4( I  5 , 3X , F 8. 4, 9X )  ) 

9694  FORMAT  ( 22X ,' SUMMARY  OF  POSITIVE  POTENTIAL  ENERGY  CHAN 
1GES  GREATER  THAN' , c7. 4 , 2X, • NT='  ,  12 , //  ,  •  ATOM     INITIA 
1L/FINAL  X     INITIAL/FINAL  Y     INITIAL/FINAL  Z     PE 
1CHANGE' /) 

9695  FORMAT  (  22X,  • SUMMAR Y  OF  NEGATIVE  POTENTIAL  ENERGY  CHAN 
1GES  GREATER  THAN  •  ,  F7  .  4, 2X  ,  'NT=  •  , I  2, //, '  ATOM     INI TI 
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1AL/FINAL    X  INITIAL/FINAL    Y  INIT I AL/ F  INAL    Z  PE 

1    CHANGE'/) 

9696  FORMAT     { 15 ,7X, 3 ( F6 . 2, 1 3X } , / , 12X, 2( F6. 2, 13X ) , F6. 2 , 
113X,F8.4) 

9697  FORMAT     ( 10X ,' INITIAL    SUMMARY    OF    POTENTIAL    ENERGY    CHANG 
1ES    ',//,('     ATOM       PE    CHANGE  ')) 

9699  FORMAT     (40X, 'FORCES    AFTER    EACH    STEP',//,'     STEP',10X, 
I'DT/STEP     • , 11X, ' FMAX' ,9X, 'ATOM     WITH    MAX    FORCE', 5X, 
l'FORCE     (1)'/) 

9700  FORMAT     ( 2Xt I  3, 6X ,  E  12.4 ,6X , E 12 . 4 , 13X , I  3 , 12 X , E12 . 4 ) 
9704     FORMAT(20Xt ' DT     CHECK    FOR    EACH    STEP ' , // t  ' STEP ', 1 1 X , 

1« STEP •*11X,»DTE1I,11X,,DTF1«,11X,«DTE,,11X, 'DTP1,//) 
9704    FORMAT( 20X, 'DT    CHECK    FOR    EACH    STE P ',//,' STE P '  ,1 1 X , 
1 'DTE1' ,11X, 'DTF1', 11X, '  DT  E  «  , 11X, 'DTF',//) 

s^sjojc^sjfsicj^^  A  3$:  A  ^c  "JzOx  ^a^csjiijc^slc^c^^^c  s^c  afc  ?*r:  3tc  >!<  a^c  afc  2^;  3^C  2^  4c 

C  'FORCE  CONFIGURATION' 

s}:  s{£  a)c  jjc  ^c  ajt  sk  ^c  #  sjc  sfc  3$:  *  #  *  #  3$c  :£  X;  *  *  *  ^  *  #  *  £  *  *  *  A  *  =fc  *  * 

c 

C      CHANGE  STATEMENT  NUMBER  9704  TO  READ 
C 

9704  FORMAT ( 20X, • DT  CHECK  FOR  EACH  STEP «,//' STEP ',  1 IX, 
l'DTF' , 12X, 'DTI  A'  ,11X, 'DTE'  ,11X,'DTIV',//) 

sjc  sfc  sfc  sjc  #  sk  #  ^c  :£  sjc  sjc  sfc  £  &  ^:  *  ik  ■%  sjc  #  A  #  5}c  *  4:  ^  *  *  #  *  *  #  *  #  *  *  ^t  sfc  s}«  sjc  #  ;£  s}<;  #  s}c  ^t  ^t  #  #  Xc  *  #  5ft  £  :£  sfc  sjc  s$c  *  ^ 

9705  FORMAT* I  5, 4( 3X , E 12 . 4)  ,/) 

C 

C      INITIALIZE  APPROPRIATE  VARIABLES  AND  VARIABLE  ARRAYS. 
C 

START=0.01*ITIME(XX) 
DO  2  I =1 ,1000 
RXK( I) =0.0 
RYK(  I  J  =  0.0 
RZK( 1 )=0.0 
VX( I )=0.0 
VY( I  )  =  0.0 
VZ(I ) =0. 3 
PKE( I )=0.0 
PPE( I)=0.0 
PTE( I) =3.0 
2  RZK  I  )  =  0.0 
ISHUT=1 
NRUN=0 

C 

C      READ  INPUT  DATA  COMMON  TO  ALL  DESIRED  CALCULATIONS. 
C 

READ  (  5,9010)  IH1 

READ  (  5,9023)  I H2 , DCON, AL PH A, RE , ROEC, ROEL 
READ  (  5,9030)  BULLE T ,SMA S ,PE XA , PEXB , I HB , THE RM 
READ  (  5,9030)  TARGET , TMAS , EXA , EXB , I HT, TEMP 
READ  (  5,9350)  IHS , PL ANE , LS , IX , I Y, I Z, CVR, MCRO  , DT I 
READ  {5,9042)  PPEPCK,PPENCK 

C 

C  INITIALIZE    CONSTANTS    AND    DEFINE     SCALING    FACTORS. 

C 
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DTIS=DTI 

R0E2=4. 0 

R0E=S0RT(R0F2) 

ROEM    =     ROE-DTI 

R0EL2  =  R0EL=*R0EL 

CVE=1.60E-19 

CVM=1.672E-27 

VFAC=0.5 

FM=1.0E-10 

FM2=FM*FM 

CVD=CVR*1. OE-10 

CVED=CVE/CVD 

PTMAS=TMAS*CVM 

PBMAS=BMAS*CVM 

HTMAS=0.5*PTr:AS/CVE 

HBMAS=0 .5*P8MAS/CVE 

TSAVE=BMAS/(BMAS+TMAS) 

BSAVE=TMAS/(BMAS+TMAS) 


C 

c 
c 


DEFINE    REPULSIVE    POTENTIAL    PARAMETERS. 


FXA=ALOG(-EXB*CVED)+EXA 
PFXA=AL0G(-PEX3*CVED)+PEXA 
PPTC=EXP(PEXA+PEXB*ROE) 
PAC=ALOG(CVED)+PEXA 

PFPTC=EXP(PAC+PEXB*ROE) 


C 

c 
c 


DEFINE  ATTRACTIVE  POTENTIAL  PARAMETERS, 


CGDl  =  ALOG(DCnN)+2.0*JH.PHA*RE 

CGD2=ALOG(2.0*DCCN)+ALPHA*RE 

CGB1=-2.0*ALPHA*CVR 

CGB2=-ALPHA*CVR 

CGF1=AL0G(-CGB1-CVED)+CGD1 

CGF2=ALOG(-CGB2*CVED)+CGD2 


C 

c 
c 


DEFINE    REGIONS     OF    APPLICABILITY    OF    POTENTIAL    FUNCTIONS 


ROEA=1.53/CVR 
ROEB=2. O/CVR 
ROEC2=ROEC*ROEC 


C 
C 

c 
c 
c 
c 
c 


DEFINE  PARAMETERS  USED  TO  DETERMINE  THE  BEST  CUBIC 
FIT  BETWEEN  THE  MAXIMUM  DISTANCE  OF  APPLICABILITY  OF 
OF  THE  REPULSIVE  POTENTIAL  AND  THE  MINIMUM  DISTANCE 
OF  APPLICABILITY  OF  THE  ATTRACTIVE  POTENTIAL. 
SUBROUTINE  CROSYM  PERFORMS  THE  NECESSARY  CALCULATIONS 


AC  1,  1)  =  1 .0 

A(1,2)=R0EA 

A(  1,3)  =ROEA=f  ROEA 

A( 1,4)=RCEA**3 

A(l»5) =EXP( EXA+EXB*ROEA) 

A(2, 1)=1.0 
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A(2,2)=R0EB 

A(2,3) =  ROEB*ROEB 

A(2,4)=R0EB**3 

A(2,5)=EXP(CGD1+CGB1*R0EB)-EXP(CGD2+CGB2*R0EB) 

A  ( 3  ♦  1 )  =  0  .  0 

A(3t 2)=-1.0 

A(3,3)=-2.0*R0EA 

A (3,4)  =-3. 0* ROE A* ROE A 

A (3,  5)=EXP(FXA+EX6*ROEA) /CVED 

A(4,l )=0.0 

A  (  4  ,  2 )  =- 1 .  0 

A(4,3)=-2.0*R0EB 

A  (4,  4)  =-3.0*R0E3*R0EB 

A(4,5) =( EXP( CG F 1+CGB1*R0EB) -EXP (CGF2+CGB2* ROE B))/ CVED 

CALL    CROSYM 

CPO=A(l ,5) 

CP1=A(2, 5) 

CP2=A(3, 5) 

CP3=A(4,5) 

CF0=-CP1*CVED 

CF1=-2.0*CP2*CVED 

CF2=-3.0*CP3*CVED 

C 

C      READ  INPUT  DATA  FOR  EACH  SITE  TO  BE  INVESTIGATED. 

C      MULTIPLE  INVESTIGATIONS  ARE  POSSIBLE  BY  SIMPLY  READING 

C       IN  DATA  FOR  MORE  THAN  ONE  SITE. 

C 

5       READ    {     5,9043)     FVR , NTT , NS , ND , I P, I  DEEP, I TY PE,     D1X, 
1D1Y»D1Z* IVACX,I  VACY,I  VACZ 
RE AD (5, 9041 , END=99  99)     DDT  I  I , DDT  I F , XNVAC , YN VAC , ZN VAC 

***  **  ************  *****  ***  *********************************** 
******&*****  **********  **  *  ***  ****  *** 

C  'FORCE    CONFIGURATION' 

*********************************** 

c 

C  IN    THE    FORCE    CONFIGURATION    ONLY,     READ    IN    DTI     VALUES 

C  TO    BE    ASSIGNED    AFTER    FORCE    COMPARISON,     AND 

C  DEFINE    VARIABLE    VFAC2. 
C 

READ (5,9  052)    DT I Al , DTI A2 ,DTI A3, DTI  VI ,DT I V2 ,DTI V3 
VFAC2=VFAC*VFAC 

*  *******************  *************************************** fc 

IF(NTT.EO.O)    GO    TO    9999 

IQ=ITYPE-1 

EV=EVR*1.0E+3 

DTI=DTIS 

TPKE=EV 

********************  *********************************X!  ****** 
*********************************** 

C  'POTENTIAL  WELL  CONFIGURATION' 

*********************************** 

C 

C  CONSTRUCT    A    'RELAXED'     CRYSTAL     IN    THE    COMPUTER    BY 

C  READING     IN     'RELAXED'    CRYSTAL    PARAMETERS     AND    POSITIONS 
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C      VICE  USING  SUBROUTINE  BlOO  TO  CONSTRUCT  THE  CRYSTAL. 
C 

READ   (5,9690)  LL , D1X, D1Y, D 1Z, NVAC 
DO  15  I =1,LL,3 
K=I  +  1 
J=I+2 
READ   (5,9691)  RX ( I )  ,  RY ( I )  , RZ ( I )  , RX ( K) , RY ( K  )  , RZ ( K) , 
1RX( J),RY(J  ),RZ I  J) 
15  CONTINUE 

C  'OTHER  CONFIGURATIONS' 

C 

C      SELECT  THE  DESIRED  CRYSTAL  STRUCTURE  AND  ORIENTATION. 

C      SUBROUTINE  BlOO  CONSTRUCTS  THE  (100)  PLANES  OF  A  BODY- 

C      CENTERED  CUBIC  CRYSTAL  IN  THE  COMPUTER. 

C 

14  CALL  BlOO 

30     ILAY=IDEEP 

IF(IDEEP)     35,35,40 
35    LD=LL 

1LAY=IY 
40    RLL=1.0/LL 

TP0TL=1 .0 

DO    45     I =1,LL 

RXK( I)=RX( I ) 

RYK( I )  =  RY( I  ) 

RZK(I ) =RZ(I ) 

RXI  (  I  )  =  RX(  I  ) 

RYK  I  )  =  RY(  I  ) 
45    RZI (I) =RZ(I) 
C 

IF(NRUN.EO.O)     GO    TO    60 

DO    55    I =1,LL 

LCUT(I )=0 

RX(  I  )  =  PXI(  I  ) 

RYU  )=RYI  (  I  ) 

RZ(I )=RZI(  I  ) 

RXK(  I)  =  RXI  (  I  ) 

RYK(  I  )  =RYI  (  I  ) 
55  RZK(  I)  =  RZI ( I  ) 
60    NRUN=1 

C 

C  THIS    SECTION    CALCULATES    THE    ENERGIES    OF    ALL    ATOMS     IN 

C  THEIR     INITIAL    POSITIONS    IN    THE     PERFECT    LATTICE     (THAT 

C  IS,     WITH    NO     INTERSTITIAL     IMPLANTED).     INITIAL    POSITIONS 

C  AND    ENERGIES    OF     ALL    ATOMS    ARE    PRINTED    TO    PRGVIDE     A 

C  COMPARISON    WITH    CHANGES    IN    CRYSTAL    ATOM    POSITIONS 

C  AND    ENERGIES    CAUSED    BY     IMPLANTATION    OF    THE     INTERSTI- 

C  TIAL. 

C 

jkskiikskskskakjjojrjk  A^  ^  ^  ^^*  *>^  ^  ^  ^*^^^*^^^^  ^^^^^^  ^sjc**^^*^*^***^:*^^*  sjc^sjcsjc* 
akstsksksfcsk*  jkskjfjkjk  skikjkjkjkskakikjkjkjk  *  skak  jkjfcsksk  aksk  ^*t 

C  'POTENTIAL  WELL  CONFIGURATION' 

*  sk  jfc  #  A  *  >k  *  *  jfc  J?  ic  *  jk  *  sk  ak  sk  *  sk  jk  >k  H  -.  jfc  si;  jfc  sfc  ;k  jp  *  *  ak  ak  ak  ak 


70 


c 

C      THIS  CALCULATION  CAN  NOT  BE  MADE  IN  THE  POTENTIAL 

C      WELL  CONFIGURATION  WITHOUT  DESTROYING  THE  INPUT 

C      DATA  OF  THE  'RELAXED'  CRYSTAL.  CONSEQUENTLY, 

C 

C  RX(1J=25.0 

C 

C      SHOULD  BE  DELETED  IN  THE  POTENTIAL  WELL  CONFIGURATION. 

C 

A************************************  *********  ************** 

RX(1)=25.0 
DO  63  1=1, LL 
VX( I )=0.0 
VYU  )=0.3 
VZ( I J=0.0 
PPE(I)=0.0 
PKE( I ) =0.0 
63  PTE(I)=0.0 
TPOT=0.0 
NEW  =  0 

CALL    ENERGY 
NPAGE=1 
NT  =  0 

WRITE     (     6,9610) 

WRITE     (     6,962J)     IH2,NT 
DO    61    1=1, LL, 3 
K=I+1 
J  =  I  +  2 

61  WRITE     (     6,9630)     I  , RX (  I )  , RY ( I  )  ,  RZ { I )  , PPE ( I  ) , K  ,RX ( K)  , 
1RY(K),RZ(K),PPE(K),J,RX(J),RY(J),RZ(J),PPE(J) 

WRITE     (     6,9650)     NPAGE 
DO    62    1=1, LL, 3 
K=I+1 
J  =  I  +  2 

PPESAV(  I  )  =  PPE(  I ) 
PPESAV(K)=PPE(K) 
PPESAVt  J)=PPE( J) 

62  CONTINUE 

*********************************************  ************:*** 
***************  ******************** 

C  'POTENTIAL  WELL  CONFIGURATION' 

*******  sje^arr;  ****rftj5c*****>!j************* 

c 

C  THE  VARIABLES  BELOW  ARE  USED  TO  CREATE  OFFSETS  FROM 

C  EQUILIBRIUM  POSITIONS  IN  THE  'RELAXED'  CRYSTAL.  IF 

C  NO  OFFSET  IS  DESIRED,  THESE  VARIABLES  SHOULD  BE 

C  INCLUDED  BUT  SET  EQUAL  TO  ZERO.  THE  OFFSET  IS 

C  INSTITUTED  BY  SUBROUTINE  PLACE.  THESE  VARIABLES 

C  SHOULD  NOT  BE  INCLUDED  IN  OTHER  CONFIGURATIONS. 
C 

D1X=0.0 
D1Y=D. 3 
D1Z=0.0 

***  **<<  ***£** ****************************************** ****** 

c 

C  SUBROUTINE    PLACE    CREATES    THE    DESIRED   VACANCY,     INTER- 

C  STITIAL,    OR     SELF     INTERSTITIAL    IN    THE    CRYSTAL. 

C  IN   THE    POTENTIAL    WELL    CONFIGURATION,     THE     INTERSTITIAL 
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C      IS  OFFSET  FROM  ITS  EQUILIBRIUM  POSITION. 
C 

CALL    PLACE 
DO    6  5    1=1, LL 
VX( I )=0.0 
VY(I )=0.0 
VZ(  I  )  =  0.0 
PPE( I  )=0.0 
PKE(  I  J  =0.0 
65    PTE(IJ=0.0 
TPOT=0.0 
NEW  =  0 

C 

C  THE    ENERGY    SUBROUTINE    NOW    CALCULATES    ENERGIES    OF    ALL 

C  ATOMS  OF  THE  CRYSTAL  AFTER  IMPLANTATION.  THESE 

C  ENERGIES  AND  THE  INITIAL  POSITIONS  OF  ALL  ATOMS  ARE 

C  PRINTED  FOR  TIME  ZERO.  CHANGES  IN  POTENTIAL  ENERGY  OF 

C  ALL  ATOMS  AS  A  RESULT  OF  IMPLANTATION  ARE  ALSO  CALCU- 

C  LATED    AND    PRINTED. 

C 


CALL    ENERGY 

BIND=-TPOT 

TE=TPOT+BIND 

TIME=0.0 

NT=0 

WRITE     (     6.962  0) 

IH2,NT 

DO    70     1=1, LL, 3 

K=I  +  1 

J  =  I  +  2 

WRITE     (     6,9630) 

I,RX( I ) 

1RY(K),RZ(K),PPE(K),J,RX( 

,RY(  I )  ,RZ(  I )  ,PPE(I  )  ,K,RX( K)  , 
J  )  ,  RY  ( J  ) ,  R  Z  (  J  )  ,  P  P  E  (  J  ) 
WRITE  (  6,9640)  TPKE , TPOT ,TE , ROE L 
NPAGE=NPAGE+1 

WRITE  (  6,965D)  NPAGE 
WRITE  (6,9697) 
DO  80  1=1, LL,^ 
K  =  I  +  1 
J=I  +  2 
L=I+3 

PPEINT1 I )=PPE( I  J-PPESAVU  ) 
PPEINTCK J=PPE(K)-PPESAV(K) 
PPEINT(J)=PPE(J)-PPESAV(J) 
PPEINT(L)=PPE(L)-PPESAV(L) 
80  WRITE  (6,9693)  I  , PP E  I  NT ( I ) , K , PPE I  NT ( K )  , J , PPEI NT ( J ) , L , 
1 PPE INT (L) 
NPAGE=NPAGE+1 
WRITE  (6,9650)  NPAGE 
DT=1.0E-15 

£#£  ji;  a}c  #  at  ;=c  :£  aj?  #  a£  afcr  £  air  a}c  ak  a*  ak  £  ate  a$r  afc  a*  ^^^^^^^^^^^^^^^^^^^^c^^ajc^c^cak^^^c^;^^^:^^* 

at  aic  afe  ate  at  >k  a&  a£c  &  a&*  ak  at*  afcakak&&&ak&?&a&a^a&afe*.&a^a&aka&2»akaka&ak 

C  "FORCE    CONFIGURATION' 

C 

C  SINCE    TIMESTEPS     IN    THIS    CONFIGURATION    ARE    DETERMINED 

C  DIRECTLY    BY    A    FORCE    COMPARISON,     NO    TIMESTEP 

C  DECREMENTING    PROCEDURE     IS    REQUIRED.    THE    ONLY   CARD 

C  NEEDED     IS 

100    CONTINUE 
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A  A  *  *  A  A  A  A  A  A  A  *  A  A  *  A  A  A  A  *  A  *  *  A  A  A  A  *  A  A  *  *  *  A  A 

C  'OTHER    CONFIGURATIONS' 

c 

C      THE  FOLLOWING  STATEMENTS  DEFINE  THE  DTI 

C      DECREMENTING  PROCEDURE. 

C 

DDTI=DDTI I 

NDEC=0 
100    DTI=DTI-DDTI 

IF(DTI.LT.DDTIF)     DTI=DDTIF 
105    NDEC=NDEC+1 

feAA'itfcAAA***********^  £  **  **  AAA*  AAA*  A*  A  A  A*  A*  A** A* A****** A  AAA*  A  A 

c 

C  THE    MAIN    BODY    OF    THE    PROGRAM    NOW    SOLVES    THE    EQUATIONS 

C  OF    MOTION    BY    THE    AVERAGE    FORCE     METHOD    AND    DETERMINES 

C  POSITIONS    0^    ALL    ATOMS     AT    THE    END    OF    THIS    TIMESTEP. 

C  SUBROUTINE     STEP    PERFORMS    ALL     FORCE    CALCULATIONS. 
C 

DTOD=DT/CVD 

TFAC=2.0*PTMAS#DTI*CVD 
TFAC8-2.0*PBMAS*DTI*CVD 
TEFAC=DTI*CVD 

HDTOD=0.5*DT0D 

DTOM=DT/PTMAS 

HDT0M=0.5*DTCM 

DTOMB=DT/PBMAS 

HDT0MB=0.5*DT0MB 
200    CALL    STEP 

IF  (LCUT (1) .GT.O)     GO    TO    240 

1  =  1 

RXK( I }=RX( I ) 

RYK(I ) =RY( I) 

RZK( I)=RZ( I ) 

RX(  I  )  =  RX{ I  )+DTOD*(HDTOMB*FX{  I  )  +  VX( I  J ) 

RY(I)=RY(I  )+DTCD*(HDTOME*FY(I  )+VY(  I  )  ) 

RZ(  I  )=RZ(  I  )+DTOD*:(HQTOMB*FZ(  I  )  +  VZ(  I  )  ) 
240     DO    245     I=2,LD 

IF(LCUT( I) .GT.O)GO    TO    245 

RXK(  I)  =  RX(I  ) 

RYKU ) =RY(I) 

RZK{  I  )  =RZ(I.) 

RX( I )=PX( I )+DTOD*(HDTOM*FX( I )+VX( I ) J 

RY(I }=RY(I )+DTGD*(HDTOM*FY( I )+VY( I ) ) 

RZ( I ) =PZ(I )+DTOD*(HDTOM*FZ(I )+VZ( I) J 
245    CONTINUE 

CALL    STEP 

FSTACC(NT)=SORT(FX( 1)**2+FY( 1 ) **2+FZ ( 1 ) **2 ) / BMAS 

EMAX=0.0 

FMAX=0. D 

TIME=TIME+DT 

NT=NT+1 

IF(LCUTtl) .GT.O)    GO    TO    265 

1=1 

VSS=VX( I ) 

VX{ I )=VSS+HDTOMB*FX(I ) 

RXt I)=RXK( I )+( VX(I )+VSS)*HDTOD 

VSS=VY(I ) 

VY( I )=VSS+HDTOMB*FY(I J 

RY(  I)  =  PYK( I )+(VY(I  )+VSS)*HDTOD 

VSS=VZ(I ) 
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VZ  (  I  )  =VSS+HDTOMB*FZ  (I  ) 

RZ(  I  )  =  PZK( I  )  +  { VZ (I  )+VSS)*HDTOD 

PKE(I)=VX(I)*VX(I )+VY(  I  )*VY(  I  )  +VZ ( I ) *VZ (  I ) 

EMAX1=PKE( I ) 

FMAX12=FX( I)*FX(I)+FY(I)*FY(I)+FZ(I J*FZ( I) 

FMAX1=SQRT(FMAX12) 

AMAX1=FMAX1/BMAS 

FSSACC(NT)=AMAXL 
260    FXCI )=0.0 

FY(  I  )  =  0.0 

FZ(  I  )  =  0.0 

FMAX=0. 0 

EMAX=0.0 

F2M=0.0 
265    DO    2  80    I =2, LD 

IF(LCUT( I ) .GT. 0)G0    TO    280 

VSS=VX( I ) 

VX( I )=VSS+HDTOM*FX(I) 

RX( I)=RXK( I )+(VX( I )+VSS)*HDTOD 

VSS=VY( I ) 

VY( I)=VSS+HDTOM*FY( I ) 

RY(  I)  =  RYK{ I )+{VY(  I  )+VSS)*HDTCD 

VSS=VZ( I } 

VZ( I )=VSS+HDTOM*FZ( I ) 

RZ( I)=RZK(I )+{VZ( I ) +VSS  )-HDTOD 

PKE{I}=VX(I)*VX(I)+VYm*VY(I)  +  VZ(I)*VZ(  I) 
275    F2=FX( I )*FX( I)+FY(I)*FY(I)+FZ(I)*FZ(I) 

F  X  {  I  )  =  J  .  0 

FY(I )  =  0. 0 

FZ(  I  )=0.0 

IF( F2.LE.F2M)    GO    TO    278 

F2M=F2 

IFMXAM{NT)=I 
278     IF(PKE( I  J .GT.EMAX)     EMAX=PKE(I) 
280    CONTINUE 

FMAX=SQRT(F2M) 

AMAXL=FMAX/TMAS 

FORMAX<NT)=AMAXL 

C  'FORCE    CONFIGURATION' 

C 

C  TIMESTEP    DETERMINATION     IN    THE    FORCE    CONFIGURATION 

C  IS    PERFORMED    BY    COMPARING    THE    MAXIMUM    FORCE     IN    THE 

C  CRYSTAL    WITH    APPROPRIATELY    CHOSEN    TEST    VALUES. 

C 

AMAX=AMAX1 

IF(AMAXL.GT.AMAX)     GO    TO    282 

IFMXAM(NT)=1 

GO    TO    284 
282    AMAX=AMAXL 
284    DTL  =  DT 

EMAXL=EMAX*VFAC2 

IF(AMAX.LE.  1.0E-8)     DTIA  =  DTIA1 

IF(AMAX.LE.1.0E-9)     DTIA=DTIA2 

IFIAMAX.LE.1.0E-10)     DTIA=DTIA3 

IF(AMAX.LE. l.OE-11)     ISHUT=-1 

CTIME=0.01*ITIME(XX)-START 

EMAX=EMAX1 

IF(EMAXL. GT.EMAX)     EMAX=EMAXL 

IF(EMAX.LE.l .OE+6)     DTIV=DTIV1 

IFCEMAX.LE. l.OE+4)     DTIV=DTIV2 

IF(EMAX.LE.1.0E+2)     DTIV=DTIV3 

DTFCK=(DTI A*CVD*2.0)/ ( AMAX/CVM) 
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IF(DTFCK.LT.O. 0) DTFCK=-DTFCK 

DTF=S0RT (DTFCK) 

DTFSP(NT)=DTF 

DTIASPl NT) =DTI  A 

IF(EMAX .GT.0.0 )     GO    TO    290 

DTE=1.0E-5 

GO    TO    295 
290    DTE=(DTIV*CVD)/SQRT(EMAX) 
295    DTESP(NT)=DTE 

DTIVSP(NT)=DTIV 

DT=DTF 

IF(DTE. LE.DT)DT=DTE 

DTSTEP(NT)=DT 

C  'OTHER    CONFIGURATIONS' 

c 

C  FOUR    NEW    TIMESTEPS    ARE  CALCULATED    BASED    ON    MAXIMUM 

C  FORCES     AND    ENERGIES    OF  INITIAL    AND    FINAL     POSITIONS. 

C  THE    SMALLEST     IS    CHOSEN  AS    THE    NEXT    TIMESTEP    INTERVAL. 
C 

DTL=DT 

CT I ME= 0.01* I  TIME (XX) -START 

DTE1=TEFAC*SQRT( 1.0/EMAX1) 

DTF1=S0RT(TFACB/FMAX1) 

DTE=TEFAC*SQRT ( 1.0/ E MAX) 

DTF=SQRT (TFAC/ FMAX ) 

DTKNT)  =  DTE1 

DT2(NT)=DTF1 

DT3(NT) =DTE 

DT4(NT)=DTF 

IF(EMAX1  .GT.EMAX)     EMAX=EMAX1 

DT=DTE1 

IF(DT.GT.DTFl)     DT=DTF1 

IF(DT  .GT  .DTE)     DT=DTE 

IF(DT.GT.DTF)     DT=DTF 

DTSTEP(NT)=DT 

c 

C  DAMPING     IS    INTRODUCED     IN    THE     FORM    OF    A    DAMPING    FACTOR 

C  WHICH    DECREMENTS    ALL    VELOCITY    COMPONENTS. 

C 

300     IF( ISHUT.EO.-l )     GO    TO    400 
310     IF(NS-NT)     403,403,320 
320    DO    325     1=1, LL 

VX( I )=VFAC*VX(  I  ) 

VY( I )=VFAC*VY( I ) 
325    VZ( I)=VFAC*VZ(  I  ) 

C  "FORCE    CONFIGURATION' 

C 

C  SINCE  THE  TIMESTEP  HAS  ALREADY  BEEN  DETERMINED  IN  THE 

C  FORCE  METHOD,  THE  COMPUTATION  CAN  NOW  SHIFT  TO 

C  STATEMENT  NUMBER  100  AND  BEGIN  THE  DYNAMICS  FOR  THE 
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C  NEXT    TIMESTEP. 

C 

GO  TO  100 

************  ********  ***  *  ********  *** 

C  • OTHER  CONFIGURATIONS' 

*********************************** 

c 

C      SHIFT  TO  STATEMENT  NUMBER  800  FOR  CALCULATION  OF  THE 

C      DTI  DECREMENT  IF  REQUIRED  FOR  THIS  TIMESTEP. 

C 

GO  TO  800 

*  **  **  *****  ****  ******  ***************  ************************* 

C 

C      SUBROUTINE  PRINT  PRINTS  PERTINENT  GENERAL  DATA  FOR 

C      EACH  TIMESTEP  PRINTOUT. 

C 

400  CALL  PRINT 

C 

C      RELATIVE  MOTION,  VELOCITY,  AND  ENERGY  OF  EACH  ATOM 

C      ARE  PERIODICALLY  PRINTED. 

C 

410  TPOT=0. 0 

DO  450  1=1, LL 

PPE(  I  )  =  0.0 
450  PTE( I ) =0.0 

CALL  ENCRGY 

PKE( 1 )=HBMAS*PKE(1 ) 

TPKF=PKF (1 ) 

PTE( 1)=PKE( 1)+PPE( 1) 

DO  620  1=2, LL 

PKE( I) =HTMAS*PKE (I ) 

TPKE=TPKE+PKE( I ) 
620  PTE( I) =PKE(I )+PPE( I ) 

TE=TPOT+BIND 
WRITE  (  6.9660) 

DTEST  =  0.1*(RX(1)-RXIU))**2 

IF  (DTEST.GT.  0.01)  DTEST=  0.01 

IF(TPOT.LE.TPOTL)  GO  TO  700 

ERAT=TPKE*RLL 
700  DO  750  1=1, LD 

DX( I )=RX (I )-RX I ( I ) 

DY( I J=RY(I J-RYI ( I) 

DZ( I )=RZ( I J-RZI ( I ) 

IF  (DX( I )**2.GE.DTEST)  GO  TO  720 

IF  (DY( I )**2.GE.DTEST)  GO  TO  720 

IF  (DZ(  I)**2.GE.DTEST )  GO  TO  720 

GO  TO  750 
C 

720   WRITE  (  6,9670)  I , DX (  I  ) , DY( I  ) , DZ ( I ) , VX( I ) , VY( I ) , 

1VZ(I ) ,PKE( I) ,PPE(I ) ,PTE(I) 
750  CONTINUE 

WRITE  (  6,96^0)  T PKE , TPOT , T E, ROEL 

NPAGE=NPAGE+1 
WRITE  (  6,9650)  NPAGE 

TPOTL=TPOT 
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IF(NT-NTT)  760,950,950 
760  DO  780  1=1. LL 

VX(  I  )=VFAC*VX(  I) 

VY(I )=VFAC~VY(  I  ) 

VZ(  I}=VPAC*VZ(  I  ) 
780  CONTINUE 

IF(ISHUT.EQ.-l)  GO  TO  950 
790  NS=NS+ND 

C  'FORCE  CONFIGURATION' 

c 

C  NO  DECREMENTING  PROCEDURE  FOR  DTI  IS  NEEDED,  SO  THE 

C  FORCE  CONFIGURATION  SHIFTS  TO  STATEMENT  NUMBER  100 

C  AND  BEGINS  COMPUTATIONS  FOR  THE  NEXT  TIMESTEP. 
C 

830  GO  TO  100 
C  'OTHER  CONFIGURATIONS' 

c 

C  THE    DTI     ALTERATION    PROCESS    IS    BEGUN    FOR    THE    NEXT 

C  TIMESTEP 

C 

800    IF  (NDEC. EQ.10)     GO    TO    810 

GO    TO    820 
810    DQTI=0.1*DDTI 

DTI  =  DTH-DDTI 

NDEC=0 
820    GO    TO    100 

950  CONTINUE 

C 

C  FINAL  POSITIONS  (IN  RECTANGULAR  COORDINATES)  AND 

C  BINDING  ENERGIES  OF  ALL  ATOMS  ARE  PRINTED  AFTER  THE 

C  LAST  TIMESTEP.  WRITE  (7,XXXX)  STATEMENTS  ARE  INCLUDED 

C  WHEN  DATA  DECKS  CONTAINING  COMPLETE  INFORMATION  FOR 

C  THE  ENTIRE  CRYSTAL  ARE  DESIRED.  ADDI T IONALL Y,  POSITIVE 

C  AND  NEGATIVE  POTENTIAL  ENERGY  CHANGES  GREATER  THAN 

C  A  PRESET  VALUE  ARE  DETERMINED  AND  PRINTED. 

C 

955   WRITE  (  6,9620)  IH2,NT 

WRITE  (7,9690)  LL , D1X , D1Y , Dl Z, NV AC 
DO  965  1=1, LL, 3 
K=I+1 
J  =  I  +  2 

WRITE  (7,9691)  RX( I ) , RY( I )  , RZ ( I  )  ,  RX( K)  , RY ( K)  , RZ ( K ) , 
1RX( J) ,RY( J),RZ(J) 
965   WRITE  (  6,9630)  I  , RX  (  I  )  , RY ( I )  , RZ ( I )  , PP E (  I  ) , K , RX ( K ) , 
1RY(K),RZ(K),PPE(K),J,RX(J),RY(J),RZ(J),PPE(J) 
WRITE  (  6,9640)  T PKE , TPOT , T E, ROEL 
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NPAGE=NPAGE+1 
WRITE  (  6,9650)  NPAGE 

DO  970  I  =1  ,LL 

PPEPOSt I J=0.0 

PPENEG(  I  )  =  0.0 

PPKEEP( I )=PPE( I)-PPESAV(I) 

IF(PPKEEP( I ) .3T.PPEPCK)  GO  TO  968 

IF(PPKEEP{ I ) .3T.PPENCK)  GO  TO  970 

PPENEGU  )=PPKEEP(I  ) 

GO  TO  970 
968  PPEPOS(  I)=PPKEEP(I  ) 
973    CONTINUE 

WRITE  (6,9694)  PPEPCK,NT 

DO  980  1=1 ,LL 

IF(PPEPOS(I ) .LT. PPEPCK)     GO    TO    980 

WRITE    (6,9696)     I , R XI (  I )  , RYI (  I )  ,RZ I ( I  )  ,RX( I  )  , RY(  I  )  , 
1RZU  )  ,PPEPOS  (  I) 
980    CONTINUE 

NPAGE=NPAGE+1 

WRITE  (6,9650)  NPAGE 

WRITE ( 6,9695)PPENCK,NT 

DO  990  1=1, LL 

IF (PPENEG(I) .GT.PPENCK)     GO    TO   990 
WRITE  (6,9696)  I , R XI (  I ) , RYI ( I )  ,  RZI  ( I )  , RX ( I )  ,  RY ( I  ) 
1RZ( I ),PPENEG( I  ) 
993  CONTINUE 

NPAGE=NPAGC+1 

WRITE  (6,9650)  NPAGE 

WRITE  (6,9699) 

DO  995  1=1, NT 
995  WRITE  (6,9700)  I , DTSTEP ( I ) , FORMAX ( I ) , I FMXAM( I ) , 
1FSSACC(  I  ) 

NPAGE=NPAGE+1 

WRITE     (6,9653)     NPAGE 

WRITE  (6,9704) 

DO  999  1=1, NT 

ifj$c$:if.%:3p1%it1x%it  ■#*£%■*: ■%??  if  ■*?%  ^e^rsjc  jk^c*^t^c^:^c^t5pyr  A«  if^if^iic-ififi^yii^.iiif^^it^at^c^ififi:^.^:^: 

C  'FORCE    CONFIGURATION' 

jjtsiojcjje  >!;  jV£  &&ifc^if  #it  if  ^&&^ifi£  ifif  if  if  $f$cif  if  if  if  %$c&& 

c 

C      DATA  PERTINENT  TO  TIMESTEP  DETERMINATION  BY  THE 

C      FORCE  METHOD  ARE  PRINTED. 

C 

999  WRITE  (6,9705)  I ,DTFSP ( I ) , DTI  A SP ( I )  ,DTESP (  I  )  ,  DT I VSP ( I ) 

%#ifc  it  &%%&±&&iz  if*  it  &&&&%%&*.*  ik&^&  if  #*<%.&■&%. 

C  'OTHER  CONFIGURATIONS' 

afe  *k  sk  if  if  if  if  jfc  if  if  sjc  yf  if  if  if  jfc  *  )k  3{C  £  &  if  #  if  if  if  sjc  5k  if  5k  #  it  if  #  if 

C 

C  DATA    PERTINENT    TO    TIMESTEP    DETERMINATION    BY    THE 

C  ENERGY    METHOD    ARE    PRINTED. 

C 

999    WRITE    (6,9705)     I , DT 1 (  I ) , DT2(  I ) , DT3( I ) , DT4( I ) 


$##:£*#£###*#£:£  #*#*$##  #*#***# 


78 


£  £  sje  *  *  *  <c  £  £  -^  if  if  if  if  jjc  -Jf.  jfc  it  Jc  *  -Jf  £  *  rfc  *  *  if  *  *  *  if  #  *  *  # 

C  'PLUCK    CONFIGURATION" 

9jc^c^^:^^3{c  ^tsj:  izifififif%iKi£ififififityfififit#-&ifi£ifixifi?if 

c 

C  SUBROUTINE  PLUCK  CAN  BE  USED  TO   PUNCH  DATA  CARDS  FOR 

C  A  SMALLER  CRYSTAL  CENTERED  ON  THE  INTERSTITIAL  FOR 

C  USE  IN  THE  DYNAMIC  PROGRAM.  DELETE  THESE  CARDS  IF 

C  DATA  FOR  THE  SMALLER  CRYSTAL  IS  NOT  DESIRED 
C 

CALL  PLUCK 

WRITE  (7,9690)  LL , D IX , D1Y , Dl Z ,NV AC , I XNEW,I  YNEW , I ZNEW 
LL  =  II 

DO  1100  1=1, LL 
1100  WRITE  (7,9691)  I , RXNE WI ( I )  ,RYNE WI  ( I  ) , RZNEWI  ( I ) , KE EP ( I ) 

jjc^c^Jj^^jjf^^fr  fc#fc  £#*#*sjefc>!t  jfcjjejjcs;::^  A  #####:£  jj:  ififif  i^-^ifif  ififif  if  *  &&ifif*z%:ifif  jfcijesfcjjcsjcajc 

1000  IF( ISHUT  )  9999  ,5,5 
9999  STOP 
END 

C 

C  SUBROUTINE  CROYSM  SOLVES  M  SIMULTANEOUS  EQUATIONS  BY 

C  THE  METHOD  OF  CROUT  IN  ORDER  TO  FIT  THE  BEST  CUBIC 

C  EQUATION  BETWEEN  THE  REPULSIVE  AND  ATTRACTIVE  PARTS 

C  OF  THE  POTENTIAL. 
C 


C 

c 


SUBROUTINE  CROSYM 


COMMON/COMA/  A(4,5),MCR0 

M=MCRO 

N  =  M  +  1 

11=1 
10D  13=11 

SUM=ABS( A( 11,11)) 

DO  120  1=11, M 

IF(SUM-ABS(A(I  ,11) )  )  110,120,120 
110  13=1 

SUM=ABS(A(I ,11  )) 
120  CONTINUE 

I  F ( 13-1 1 )  130,150,130 
130  DO  140  J=l  ,N 

SUM  =  -A(  II,  J') 

A(  II, J )  =  A(  13, J  ) 
140    A(I3,J)=SUM 
150    13=11+1 

DO    160     1=13, M 
160    A(  I  ,I1)=A(  I ,11  )/A(Il,I 1) 
170    J2=I1-1 

13=11+1 

IF(J2)     180,200,180 
180    DO    190    J=I3,N 

DO    190    1=1, J2 
190    A(I1,J)=A(I1,J)-A(I1,I)*A(I,J) 

IF(Il-M)     200,220,200 
200    J2=I1 

II  =  1  1+1 

DO    210     1=11, M 
DO    210    J=1,J2 
210    A(I,I1)=A(I,I1)-A(I,J)*A(J,I1) 
IF(Il-M)     100,170,100 
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•  220    DO    240    I  =1,M 

J2=M-I 

I3=J2+1 

A(I3tN) =A(  I3,N)/A(  13,13) 

IF(J2)     230,250,230 
230    DO    240    J  =  l, J2 

2  40    A(J,N)=A(J,N)-A(I3,N)*A(J,I3) 
250    RETURN 

END 

C 

C  SUBROUTINE    B100    GENERATES    A    BODY-CENTERED    CUBIC 

C  LATTICE     IN    THE     (100)     DIRECTION. 

C 

c 

SUBROUTINE    B100 


COMMON/ COM 1/RX (500) ,RY(500),RZ( 500),LCUT{ 5  00) , 
1LL,LD,ITYPE,NVAC 

C0MM0N/C0M4/IX,IY,IZ,SCX, SCY,SCZ,IDEEP,D1X,D1Y, 
1D1Z,IVACX,IVACY,IVACZ 

DIMENSION    YLAX(20) 

DATA    YLAX/20*0.0/ 

YLAX(l) =-0.20 

YLAX( 2)=-0.03 

SCX=1.0 

SCY-1.0 

SCZ=1.0 

M    =    2 

JT=0 

Y=-SCY 

DO    60    J=l, IY 

Y=Y+SCY 

KT  =  0 

Z=-SCZ 

DO    59    K=1,IZ 

Z=Z+SCZ 

IT=0 

X=-SCX 

DO    58    1=1, IX 

x=x+scx 

IF(IT-(IT/2)*2)     21,11,21 

11  IF(JT-( JT/2)*2)     57,12,57 

12  IF  (KT-(KT/2)*2)     57,30,57 

21  IF( JT-( JT/2)*2)     22,57,22 

22  IF(KT-(KT/2)*2)     30,57,30 
30    RX(M)=X 

RY(M)=Y+YLAX( J) 

RZ(M)=Z 

M  =  M+1 

IF    ( IT.NE.I VACX)    GO    TO    57 

IF    ( JT  .NE.IVACY)    GO    TO    57 

IF     (KT.NE. IVACZ)     GO    TO    57 

NVAC=M-1 

57  IT=IT+1 

58  CONTINUE 
KT=KT+1 

59  CONTINUE 

JT  =  JT  +  1 
IF(IDEEP-JT)  60,110,60 

60  CONTINUE 
LL=M-1 

100  RETURN 
110  LD=M-1 

GO  TO  60 

END 
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*******:****************************************************£ 
****  *****************  *>?*  a^ajcsjcs^t  ****  *** 

C  'POTENTIAL  WELL  CONFIGURATION' 

****  ******************************* 

c 

C  SUBROUTINE    PLACE    OFFSETS    THE     INTERSTITIAL     IN    THE 

C  'RELAXED'     CRYSTAL. 

C 


SUBROUTINE  PLACE 
C 
C 

COMMON/COM 1/RX( 500 ) , R Y( 5  00) , RZ ( 500) »LCUT(5  00) t 
1LL,LD,  ITYPE,NVAC 

COMMON /C0M3/R XI ( 50  0) , RYI ( 5  00 ) , R Z I ( 500 ) ,CVR,EVR, 
1NT,TIME,DT,DTI , ILAY,RXK( 5  00) ,RYK( 500) ,RZK(500) 

C0MM0N/C0M4/IX, IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y, 
lDl'Z,  TVACXtl  VACY,I VACZ 

C0MM0N/C0M9/XNVACYNVACZNVAC 

GO    TO       ( 10 , 20, 3  0,4  J) ,     I  TYPE 
10    LCUT(NVAC)     =    1 

LCUTQ  )     =    1 

RX(1) =0.0 

RY( 1)=0.0 

RZ(1 J =0.0 

GO    TO    50 
20    RX(1)  =  RX(  D+D1X 

RY  (1)=RY(1)+D1Y 

RZ(1)=RZ(  D+D1Z 

GO    TO    50 
3D    LCUT(NVAC)     =    1 

RX(1)     =    RX(MVAC) 

RY(1)     =     RY(NVAC) 

RZ(1)     =    RZ(NVAC) 

GO    TO    50 
40    RX(1 )  =  RX(1  )  +  DlX 

RY(  1)=P  Y(  D+D1Y 

RZ( 1>=RZ(1)+D1Z 
50    CONTINUE 

RX(NVAC)  =RX(  NVAO  +  XNVAC 

RY(NVAC)=RY(Ny/  AO  +  YNVAC 

RZ  (NVAC)  =  RZ(N\/AC)  +  ZNVAC 

RXI (NVAC)=RX(NVAC) 

RYI (NVAC)=RY(NVAC) 

RZI(NVAC)=RZ(NVAC) 

RXK(NVAC)=RX(NVAC) 

RYK(NVAC)=RY(NVAC) 

RZK(NVAC)=RZ(NVAC) 

RXI (1) =RX(1) 

RYI  (  1)=RY(  1) 

RZI  (1 )  =  RZ(1 ) 

RXK(l) =RX(1) 

RYK(1)=RY(1) 

RZK(1)=RZ(1 ) 

RETURN 

END 

*********************************** 

C  'OTHER    CONFIGURATIONS' 

*********************************** 

C 

C  SUBROUTINE    PLACE    CREATES    A    VACANCY    OR    IMPLANTS    AN 
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C  INTERSTITIAL    OR     SELF    INTERSTITIAL    IN    THE     LATTICE. 

C  LATTICE     ATOMS-   ARE    ALSO    OFFSET    AS    REQUIRED    BY 

C  INPUT    DATA. 
C 


SUBROUTINE    PLACE 
C 

c 

COMMON/ C0M1/RX(  500) ,RY(500) ,RZ( 500)  ,LCUT( 500)  , 
ILLtLD, ITYPE,NVAC 

C0MM0N/C0M3/RXI  ( 50  0)  » RYI  ( 5  00 ) , R Z I ( 500 ) , C VR  ,  EVR, 
1NT,TIME,DT,DTI , I L AY , R XK ( 5  00) ,RYK( 500) ,RZK(500) 

C0MM0N/CCM4/IX, IY, I Z , SCX , SCY , S CZ , I DE EP t D IX , D 1Y, 
1D1Z, IVACX,I VACY,I VACZ 

C0MMGN/C0M9/XMVAC»YNVACtZNVAC 

GO    TO       (13,20,30.43),     ITYPE 
10    LCUT(NVAC)    =    1 

LCUTd  )     =    1 

RX( 1) =0.0 

RY(1)=0.0 

RZ(1 )=0.0 

GO    TO    50 
20    RX(1)=RX(NVAC)+D1X 

RY(1)=RY (NVAC) +  D1Y 

RZ(1)=RZ(NVAC)+D1Z 

GO    TO    50 
30    LCUT(NVAC)     =    1 

RX(1)     =    RX(NVAC) 

RY(1)     =     RY(NVAC) 

RZ(  1)     =    RZ(NVAC) 

GO    TO    50 
40     RX(1  )  =  RX  (NVAO+D1X 

RY( 1)=RY(NVAC) +D1Y 

RZ(  1)  =  PZ  (NVAO+D1Z 
50    CONTINUE 

RX(NVAC) =RX( NVACJ+XNVAC 

RY(NVAC)=RY(NVAC)+YNVAC 

RZ(NVAC)  =  RZ(W  AU  +  ZNVAC 

RXI (NVAC)=RX(NVAC) 

RYI (NVAC)=RY(NVAC) 

RZI (NVAC)=RZ(NVAC) 

RXK(NVAC)=RX(NVAC) 

RYK (NVAC )=RY (NVAC) 

RZK(NVAC)=RZ(NVAC) 

RXK  1)=RX(  1) 

RYI  (1 )=RY(1 ) 

RZK  1)=RZ(  1) 

RXK( 1)=RX( 1) 

RYK(l) =RY(1) 

RZK( 1)=RZ( 1) 

RETURN 

END 

sjcjfcik  j)c  #5$:*  #*:£**  £**#:£  s^jjc  +  ajc  *#:$:  $*  *#*  *?.%:%%:$:'&*&***  *****  *  **  *  ********** 

c 

C  SUBROUTINE    STEP    CALCULATES    FORCES    ON    THE     INTERSTITIAL 

C  AND    ALL    OTHER    ATOMS     OF    THE    CRYSTAL. 

C 

C 

SUBROUTINE  STEP 


C0MM0N/C0M1/RX( 500) ,RY(500) , RZ(500) ,LCUT( 5  00)  , 
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100 

101 

102 

105 

110 

113 
117 
120 

123 
127 
130 

133 
137 
140 

150 

160 
162 

165 

167 

180 
190 


195 
200 
205 

210 

213 
217 
220 

223 
227 
230 

233 
237 

2  40 

250 

255 
260 


LL,LD, I 

COMMON/ 

IOtTSAV 

COMMON/ 

COMMCN/ 

CF0.CF1 

IF    ( 10- 

IP=2 

GO    TO    2 

IP  =  1 

GO    TO    2 

1=1 

IP  =2 

DO    195 

IFtLCUT 

DRX=RX( 

IF(DRX) 

IF(DRX+ 

IF (DRX- 

DRY=RY( 

IF (DRY  ) 

IF(DRY+ 

IFCDRY- 

DRZ=RZ( 

IF(DRZ) 

IF(DRZ+ 

IF (DRZ- 

DIST=DR 

IFIDIST 

DIST=SO 

IF(DIST 

FORCE=E 

GO    TO    1 

DFF=ROE 

IF(DFF- 

FORCE={ 

IF(FM-F 

FOD=FOR 

FA=FOD* 

FX( J  )  =  F 

FX(  I  )=<= 

FA=FOD'" 

FY( J  )  =  F 

FY(I )=F 

FA=FOD* 

FZ( J)=F 

FZ(I )=F 

CONTINU 


TYPE, NV AC 

CCM5/R0E,R0E2,R0EM,EXA,EXB,PEXA,PEXB,FXA,PFXA, 

EtBSAVE 

C0M6/FX( 500) ,FY( 500) ,  FZ( 500) , PAC , P FPTC  ,F M 

C0M8/R0EA,R0EB,R0EC,R0EC2,  CPO , CP1 , CP2, CP3, 

, CF2,C3D1,CGD2,CGB1 ,CGB2 ,CGF1 ,CGF2 

1)    100,101,102 

00 

00 


J  =  IP 
(J)  ) 
J)-R 

113 
ROE  ) 
ROE) 
J)-R 

123 
ROE) 
ROE) 
J  )-R 

133 
ROE) 
ROE) 
X*DR 
-ROE 
RT(D 
-ROE 
XP(P 
80 

-DIS 
1  .OE 
EXP( 
CRCE 
CE/D 
DRX 
X(J) 
X(  I  ) 
DRY 
Y(J  ) 
Yd  ) 
DRZ 
Z(J  ) 
Z(I  ) 
E 


,LL 
195,  110,195 

X(I  ) 

, 117,  117 
195, 195,120 
120, 195,195 

Yd  ) 

,127, 


127 
195 

195 


195,195,133 

130,  195,195 
Z(I  ) 
,137,137 

195, 195,140 

140,  195,195 
X+DRY~DRY+DRZ*DRZ 
2)     150,195,195 
1ST) 

M)     16  2,162,165 
FXA+PEXB*DIST) 

T 

-10)     195,195,167 

PAS+PEXB*DI ST)-PFPTC)/DFF 

)     190,190,195 

1ST 

+  FA 
-FA 

+  FA 
-FA 

+  FA 
-FA 


DO  300  I  =1 
IF(LCUT( I) 
IP=I+1 
DO  295  J=I 
IF(LCUT( J) 
DRX=RX( J)- 
IF(DRX)  21 
IFIDRX+ROE 
IF (DRX-ROE 
DRY=RY( J )- 
IF (DRY)  22 
IF  (DRY+ROE 
IF(DRY-ROE 
DRZ=RZ( J )- 
IF(DRZ)  23 
IF(DRZ+ROE 
IF (DRZ-ROE 
DIST  =  DRX*D 
IF(DIST-RO 
DIST=SORT( 
IF(DIST-RO 
IF(DIST-RO 
FORCE=EXP( 


P,LD 

)    300,205,300 


P,LL 

)    29 

RXd 

3,21 

C)     2 

C)    2 

RY(I 

3,22 

C)    2 

C)    2 

RZd 

3,23 

C)    2 

C)    2 

RX+D 

EC2) 

DIST 

EA) 

EB) 

FXA* 


5,21 
) 

7,21 

95,2 
20,2 
) 

7,22 
95,2 
30,2 
) 

7,23 
95,2 
40,2 
RY*D 
250 
) 

260, 
265, 
EXB* 


0,295 

7 

95,220 

95,295 


95,230 
95,295 


95,240 
95,295 
RY+DRZ*DRZ 
,295,295 

255,255 
270, 270 
DIST) 
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GO    TO    280 
2  65    F0RCE=DIST*(DIST*CF2+CF1)+CF0 

GO    TO    280 
2  70    F0RCE=EXP(CGF1+CGB1*DI ST )-E XP ( CGF 2+CGB2*D I  ST) 
280    F0D=F0RCE/DIST 

FA=FOD*DRX 

FX( J)=FX(J )+FA 

FX(I  )=FX(I )-FA 

FA=FOD*DRY 

FY(J  )=FY(J)+FA 

FY(I )=FY(I )-FA 

FA=FOD=?-DRZ 

FZ( J  )  =  FZ( J )+FA 

FZ(  I  )  =  FZ(  I  )-FA 
295    CONTINUE 
300    CONTINUE 

RETURN 

END 


C 

c 
c 
c 
c 


SUBROUTINE  ENERGY  CALCULATES  MUTUAL  POTENTIAL  ENERGIES 
FOR  THE  INTERSTITIAL  AND  ALL  OTHER  ATOMS  OF  THE 
CRYSTAL. 


C 

C 


SUBROUTINE    ENERGY 


COMMON /COM 
1LL,LD, ITYP 

COMMON/COM 
1IQ,TSAVE, B 

CCMMON/CCM 

CCMMON/CCM 
1CF0,CF1,CF 

IF     (10-1) 
100    IP=2 

GO    TO    200 
1D1    IP=1 

GO    TO    200 
102    1=1 

IP=2 
105    DO    595    J=I 

IF (LCUT (J) 
510    DRX=RX(J)- 

IF(DRX)  51 
513  IF(DRX+ROE 
517  IF(DRX-ROE 
520    DRY=RY(J)- 

IF(DRY)  52 
523  IF(DRY+ROE 
527  IF (DRY-ROE 
530    DRZ=RZ(J)- 

IF(DRZ)  53 
533  IF(DRZ+ROE 
537  IF(DRZ-ROE 
540    DIST=DRX*D 

IFIDIST-RO 
550  DIST=SORT{ 
560  POT=EXP(PE 
580    TPOT=TPOT+ 

PPE( I )=PPE 

PPE (J) =PPE 
595  CONTINUE 
600  CONTINUE 


1/RX( 500) ,RY(500) »RZ( 500) ,LCUT(500) » 

C,N\/ AC 

5/R0E,R0E2tR0EM,EXA,EXB,PEXA, P EXB , FXA, PFX At 

SAVE 

7/PPTC,TPOT,PPE(1000J » T LP E , ROEL , R0EL2, NEW 

8/ROEA,ROEB,ROEC,ROEC2tCP0.CPl,CP2,CP3, 

2.CGDl,CGD2tCGBl,CGB2tCGFl,CGF2 

100,101,102 


P»LL 
)     595, 
R  X  ( I  ) 
3,517, 
)     595  , 
)     52  0, 
RY(I  ) 
3,527, 
)     595, 
)     530, 
RZ(I  ) 
3,537, 
)     595  , 
)     540, 
RX+DRY 
E2)     55 
DIST  ) 
XA+PEX 
POT 

(  I  )  +  BS 
(JJ+TS 


510,595 

517 

595,523 

595,595 

527 

595,530 

595,595 

537 

595  ,540 

595,595 

*DRY+DRZ*DRZ 

0,595,595 

B*DIST  )-PPTC 

AVE*POT 
AVE*POT 


200    DO    300    I=IP,LD 
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205 


210 

213 
217 
220 

223 
227 
230 

233 
237 

2  40 

250 

255 

260 

265 

270 
280 


295 
300 


IF(LCUT 

IP=I+1 

DO    295 

IF(LCUT 

DRX=RX( 

IF(DRX) 

IFIDRX+ 

IF(DRX- 

DRY=RY( 

IF (DRY  ) 

IF(DRY+ 

IFCDRY- 

DRZ=RZ( 

IF(DRZ) 

IFCDRZ+ 

IF1DRZ- 

DIST=DR 

IF(DIST 

DIST=SO 

IFCDIST 

IF (DIST 

POT=EXP 

GO    TO    2 

POT=DI S 

GO    TO    2 

POT=EXP 

TPOT=TP 

SAVE=0. 

PPE(  I)  = 

PPE( J) = 

CONT  INU 

CONTINU 

RETURN 

END 


(I))    300,205,300 


J  =  IP 
(J)  ) 
JJ-R 

213 

ROEC 
ROEC 
J  )-R 

223 
ROEC 
ROEC 
JJ-R 

233 
ROEC 
ROEC 
X*DR 
-ROE 
RT(D 
-ROE 
-ROE 
(EXA 
80 

T*(D 
80 
(CGD 
OT+P 
5*P0 
PPE( 
PPE( 

E 

c 


,LL 

295 
X(I  ) 
,217 
)  29 
)  22 
Yd  ) 
,227 
)  29 
)  23 
Z(I  ) 
,237 
)  29 
)  24 
X+DR 
C2) 
1ST) 

A)  2 

B)  2 
+  EXB 


,210,295 

,217 

5,295,220 
0,295,295 

,227 

5,295,230 

0,295,295 

,237 

5,295,240 
J, 295, 295 
Y*DRY+DRZ*DRZ 
250,295,295 

60,255,255 
65,270, 270 
*DIST) 


IST*(DI ST*CP3+CP2)+CP1)+CP0 

1+CGB1*DIST)-EXP(CGD2+CGB2*DIST) 

OT 

T 

I )+SAVE 

JJ+SAVE 


C 
C 

c 

c 


SUBROUTINE    PRINT     PRINTS 
EACH    TIMESTEP    PRINTOUT. 


PERTINENT    GENERAL    DATA    FOR 


C 
C 


SUBROUTINE  PRINT 


COMMON 
1LL.LD, 

COMMON 
1TARGET 
1DDTIF 

COMMON 
INT, TIM 

COMMON 
1D1Z, IV 

COMMON 
110, TSA 

COMMON 

1CF0,CF 

9710    FORMAT 

9720    FORMAT 

1    UNIT 
9730    FORMAT 
IE    TEMP 
1V/) 

9740  FORMAT 
1  F5.2, 
1     )  ,  ,    4 

9741  FORMAT 
1    F5.3, 


/C0M1/RX 
ITYPE,NV 
/C0M2/H 
(4),TMAS 

/C0M3/RX 
E,DT,DTI 
/C0M4/IX 
ACX,IVAC 
/C0M5/R0 
VE,BSAVE 
/C0M8/R0 
1,CF2,CG 
(  4  OX  ,  1  OA 
(9H  TAR3 
=,F7.4,4 
( 4X, 6HMA 
=F5.2,7 

(2H  (,A4 
21HKEV, 
X,  16HVA 
(2H  { ,A4 
21HKEV, 


( 500) ,RY(500) ,RZ( 500) ,LCUT( 5  00) , 

AC 

1(20)  ,IH2(8)  ,IHS(10) , IHB(6) 

, BULLETt 4)  ,BMAS, PLANE, TEMP  , 


, IHT(6) , 
THERM, DDTI I, 


I(50  0),RYI(5  00),RZI(500),CV 
, ILAY,RXK(500) ,RYK{ 500) ,RZK 
, IY,  IZ,SCX,SCY,SCZ, IDEEP,D1 
Y, I VACZ 
E,R0E2,R0EM, EX  A, E XB , PE XA, P E 


R,EVR, 
(  500) 
X,D1Y, 

XB,FXA,PFXA, 

CP2,CP3, 


PA,ROEB,ROEC,ROEC2,CPO,CP1 , 

D1,CGD2,CG31,CGB2,CGF1,CGF2 

4,/,28X,20A4,/) 

ET  -,4A4,10HPRIMARY  -  ,4A4 , 1 X , 14HLATT I CE 

H  ANG) 

SS  =,F7.2,13X,6HMASS  =,F7.2 

H  DEG  K,,18H   THERMAL  CUTOF 


,9X,14HLATTIC 
F  =,F5.2,3H  E 


,8H)  PLANE, ,18H 

CRYSTAL  SIZE  ( 

CANCY  IN  SITE  , 

,8H)  PLANE, ,18H 

CRYSTAL  SIZE  ( 


PRIMARY 
, 12, 3H  X 
14/) 

PRIMARY 
,  I2,3H  X 


ENERGY  =, 

t I2f 3H  X  ,  12, 3H 

ENERGY  =, 

, 12, 3H  X  , 12, 
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13H    ),,    4X,I3,'     LAYERS    ARE    FREE    TO    MOVE',/) 

9742     FORMAT (2H     (,A4,8H)     PLANE, »  18H       PRIMARY    ENERGY    =, 

1    F5.2.21HKEV,        CRYSTAL    SIZE     (     ,I2,3H    X    ,I2,3H    X    ,I2,3H 
1     ),,     4X,     20HREPLACEMENT    IN    SITE     ,     14/) 

9750     FORMAT <•      IMPLANTATION    AT    SITE    #     ' , I  3, 3X,  ' X= • , I  2 , 4X, 
l'Y  =  '  ,I2,4X,'Z  =  '  ,  I2,5X, • 10='  ,I2,4X,  «DDTI 1  =  '  ,F6.4,4X, 
l'DDTIF=« ,F6.4, //, •     LATTICE    ATOM    START    POINT', 5X, 
l'X=' ,F5.2,3X,« Y=', F5.2,3X, •  Z=  • , F5.2,6X, ' INTERSTITIAL 
1  START    POINT'  ,5X, • X='  ,F5.2,3X,'  Y  =  »  ,F5.2,3X,' Z  =  '  , 
1F5.2,.// ) 

9760    F0RMAT(12H    POTENTIAL       , 6A4, 3X, 5  HP  EX A= , F9 . 5 , 2X , 5HPEXB= , 
1F9.5.2X, 5HPFXA=,F9. 5) 

9765    FORMAT ( 12X, 6A4t 3X, 5HEXA    = , F9. 5, 2X, 5HEXB    =  ,F9. 5 ,2X ,5HFX 
1A    =  ,F9.5/) 

9770    FORMATt  ■     WHEN',F8.4,'     <    R    <»,F8.4,'     THE     MATCHING    POTEN 
1TIAL    PARAMETERS     ARE',//,*  CPO    =',F10.3,',    CP1    =• 

1F13.3,',    CP2    =',F10.3,',    CP3    ='  ,  F10. 3, / , '  CFO    =' 

1E10.3,',    CF1    =',E10.3,',    CF2    =',E10.3,//) 

9780    FORMAT!1     CUT-3FF    AT',F5.2,',     WHEN    R    >    »,F6.3,«    LU,    MOR 
1SE    POTENTIAL    PARAMETERS    ARE',     8A4,//,10X,'  CG01     =', 

1F8.4,',     CGD2    =',F8.4,',     CGB1    =',F8.4,',     CGB2    =',F8.4, 
1',    CGF1     =  ',F8.4,»,     CGF2    =',F8.4,//) 

9791     FORMAT ( 36X, »DDTII=',F6.4/,36X,«DDTIF=',F6.4) 
WRITE     (     6,9710)     IHS,IH1 
WRITE     (     6,9720)     TARGET , BULLET , CVR 
WRITE     (     6,9730)     TMA S , BMA S , TE MP , THERM 
GO    TO     (401,402,403,402),     ITYPE 

401  WRITE     (     6,9740)     PLANE , E VR , I  X, I Y , I Z , NVAC 
GO    TO    405 

402  WRITE     (6,9741)     PLAN E, EVR, IX , I Y , I Z , I L AY 
GO    TO    405 

403  WRITE     (     6,9742)     PL ANE , E VR , I  X , I Y , I Z ,NVAC 
WRITE     (     6,9763)     I  HB,  PEX  A  ,  PEX  B  ,  P  FX  A 
WRITE     (     6,9765)     IH T  ,E XA , EXB , FXA 

WRITE     (     6,9770)     R OEA, ROEB, C PO , CP  1, CP 2, CP 3 , C F 0, CF 1 , CF2 
WRITE     (     6,978D)     ROEC , ROEB , I H2 , CGD1 , CGD2 , CGB1 , CGB2 , 
1CCF1    CCF? 

WRITE     (    6,9790)     NT  ,DT I , T I  ME , DT 
RETURN 
END 


C  'PLUCK    CONFIGURATION' 

C 

C       SUBROUTINE  PLUCK  CHOOSES  THE  ATOMS  WHICH  WILL 

C  MAKE  UP  THE  SMALLER  CRYSTAL  USED  IN  INITIAL  DYNAMIC 

C  CALCULATIONS.  PERTINENT  DATA  FOR  THIS  SMALLER 

C  CRYSTAL  AND  POSITIONS  OF  ALL  ATOMS  OF  THIS  CRYSTAL 

C  CAN  THEN  BE  PUNCHED   CUT  ON  DATA  CARDS  TO  BE  USED  AS 

C  INPUT  TO  THE  DYNAMIC  PROGRAM. 


C 


SUBROUTINE    PLUCK 

C0MM0N/C0M1/RX( 500)  ,RY(500)  ,RZ(500)  ,LCUT(5  00)  , 
1LL,LD,  ITYPE,  N\/AC 

C0MM0N/COM4/IX,I Y,I Z, SCX , SCY ,S C Z, I  DEEP, Dl X , D1Y ,  Dl Z, 
1IVACX,  IVACY,I\/ACZ 

C0MM0N/C0M10/IXNEW,  I YNEW , I ZNEW ,  1 1 

C0MM0N/C0M11/RXNEWI (250),RYNEWI(2  50),RZNEWI(250), 
1KEEP(250),NNUM(250) 
1500     IXNEW=7 

IYNEW=IVACY+3 
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IZNEW=7 

NM=5 

NI  =  8 

11  =  2 

KM  =  0 

NX=16 

NII3=10 

NI 14=5 

IF(  IYNEW.E0.3)     GO    TO    1514 

IF(IYNEW.E0.5)     GO    TO    1514 
1505    DO     1539    1=1  I • NM 

NNUM(  I)=NI 
1509    NI=NI+1 

NI =NI  +  1 

I  I  = I  1+4 

NM=NM+4 

IF    (II. LE.NX)    GO    TO    1505 

NX=NX+9 

NI=NI+NI 14 

MM=MM+1 

IF( IYNEW.EQ.MM)     GO    TO    1600 

NM=NM-1 

GO    TO    1515 

1514  NX=9 
NM=4 
NI 13=4 
NII4=11 

1515  DO    1520    I=II,NM 
NNUM( I) =NI 

1520    NI=NI+1 

N I  =N I  +  2 

11=11+3 

NM=NM+3 

IF( II. LE.NX)     GO    TO    1515 

NX=NX+16 

NI=NI+NI 13 

MM=MM+1 

IF( IYNEW.EQ.MM )    GO    TO    1600 

NM=NM+1 

GO    TO    1505 
1600     11=11-1 

RXNEWI { 1)=RX(1 ) 

RYNEWH  1)  =  RY(1  ) 

RZNFWI (1  )  =  RZ(1 ) 

KEEP( 1) =1 

NNUMC 1 )=1 
1700    DO    1750    1=2,11 

RXNEWI ( I )=RX(NNUM( I ) ) 

RYNEWI (  I  )=RY(NNUM(  I  )  ) 

RZNEWI (I )=RZ(NNUM( I ) ) 

KEEP( I )=NNUM(I ) 
1750    CONTINUE 

RETURN 

END 


BLOCK    DATA 

C0MM0N/CCM1/RX( 1000) , RY(IOOO) ,RZ(1000) ,LCUT(1000) , 
1LL,LD, ITYPE,NVAC 

DATA    RX/ 100  0*0.0/, RY/1 000*0.0/ » RZ/ 1000*0 .0/ t 
1LCUT/1000*0.  0/ 

COMMON/ C0M3/RX I  ( 10  00)  , RYI (  1000 )  ,  RZ  I  (  1000)  ,CVR  ,EVR , 
1NT,TIME,DT,DTI,ILAY,RXK(13  0  3),RYK(1000),RZK(1000) 

DATA    RXI/1000*0.0/,RYI /l 000*0. 0/ .RZ I /l 000*0. 0/ 

COMMON/ CCM6/FXI 1000),  FY(  1000  ),FZ(  1000) , P AC , P FP TC , FM 

DATA    FX/103  0*0.0/,FY/1300*0.0/,FZ/1000*0.0/ 

END 
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-  COMPUTER  PROGRAM 
FOR  DYNAMIC  SIMULATIONS 


THE 
ENERGY  RED 
STITIAL  AT 
REQUIRED  A 
SIMULATION 
TO  POINT  0 
PROGRAM.  T 

(  1)  T 
ATOM  CRYST 

(  2)  T 
MINED  BY  S 


DYNAMIC  PR 
UIRED  IN  b 


OM 
RE 
S. 
JT 
WO 
HE 
AL 
HE 


TO  EXIT 
ESSENTI 
CONSEOJ 
SIGNIFI 
CONFIGJ 
ENTIRE 
IS  USED 
PLUCK  Z 


U3ROUTINE 


OGRAM    IS    USED    TO    DETERMINE    THE    MINIMUM 
SPECIFIC     DIRECTION    TO    CAUSE    AN     INTER- 
FROM    THE    CRYSTAL.     THE    COMPUTATIONS 
ALLY    THE    SAME    AS    THOSE     IN    THE    STATIC 
ENTLY, COMMENTS    WILL    ONLY    BE     INCLUDED 
CANT     DIFFERENCES    FROM    THE     STATIC 
RATIONS    OF    THE    DYNAMIC    PROGRAM    EXIST: 
CRYSTAL    CONFIGURATION    -    THE     ENTIRE    250 

IN    ALL    CALCULATIONS. 
ONFIGURATIGN    -    A    SMALLER    CRYSTAL    DETER- 
PLUCK    IS    USED    FOR    ALL    CALCULATIONS. 


D I  MENS  I 
D I  MENS  I 
DIMENS  I 
DIMFNSI 
COMMON/ 
LL,LD, I 
COMMON/ 
TARGET! 
COMMON/ 
NT, TIME 
COMMON/ 
COMMON/ 
IO,TSAV 
COMMON/ 
COMMON/ 
COMMON!/ 
CF0,CF1 
COMMON/ 
COMMON/ 


ON    VX 

ON  DX 
ON  RX 
ON  KE 
C0M1/ 
TYPE, 
COM  2/ 
4  )  ,  T  M 
C0M3/ 
,DT,D 
COM  4/ 
C0M5/ 
EtBSA 
COM  6/ 
COM  7/ 
COM  8/ 
,CF2, 
COM  9/ 
COMA/ 


(1000) 
(1000) 
K( 1000 
EP  (500 
RX( 100 
N\/AC 
UK  20 
AS,BUL 
RX  I  ( 1  0 
TI  ,ILA 
IX, IY, 
ROE.RO 
VE  ,CX, 
FX (100 
PPTCT 
ROEA,R 
CGD1,C 
RXSAVE 
A(^,5 


•VY(1303)tVZ(10l 

,DY(  1000)  ,DZ( 101 

) ,RYK(  1000), RZK 

) 

0),RY(1000),RZ(  1000)  , LCUT(IOOO)  , 


00),PKE(  1000) 

100) ,PTE(1000) 
(  1000) 


)  ,IH 

LET( 
00), 
Y 

IZ,  S 
E2,R 
CY,C 
0),  F 
POT, 


2(8) ,IHS (10 
4) ,BMAS,PLA 
RYK  100D  ),R 

CX,SCY,SCZ, 

OEM, EXA,EXB 


Y( 1000  ), FZ( 
PPE( 1000)  ,T 
R0ECR0EC2, 
CGB1,CGB2,C 
AVE,RZSAVE 
) ,MCRO 


u:di 
GC2, 
,  RYS 


)  ,IHB(6) ,  IHT(6), 

NE, TEMP, THERM 

ZI ( 1000 ) ,CVR,EVR, 

IDEEP ,D1X,D1Y,D1Z 
,PEXA, PEXB, FXA,PFXA, 

1000) ,PAC,PFPTC,FM 
LPE,ROEL, R0EL2,NEW 
CP0,CP1 ,CP2,CP3, 
GF1,CGF2 


9010 
9020 
9030 
9040 
9050 

9610 

9620 

9630 
9640 


FORMATt  20A4) 
FORMAT (8A4,3F8 
F0RMAT(4A4,3F3 
FORMAT! F6.5, 5X 
F0RMATU0A4,  A4- 


5,2F5.2) 
5,6A4,F6 .2) 
I5,6I4,3F5.2) 
413, F8  .4,  I4,F5.3) 


NT    =  ■ 14,//, 
•),//) 


9650 
9660 
] 
9670 
9680 


FORMATt  1H1) 

FORMAT (47X, 'SJMMARY  OF  ATOMS '//, 35X, 8A4, 

3(  •  ATOM       POSITION       BIND  ENERGY 

F0RMAT(3(I5,3F6.2,F8.4,8X)) 

FORMAT(/4X,F10 .3,25H  EV, TOTAL  KINETIC  EN ERGY, , F 10 . 3 , 
127H  EV, TOTAL  POTENTIAL  ENE RGY , F10. 3 , '  EV ,  REDUCTIONS/ 
1/60X, 'RADIUS  = ' , F5.2, ) 

FORMAT  (U5X,4H PAGE,  13,  /,1H1  ) 

FORMAT(/   '      ATOM     DX         DY 

VX         VY         VZ        KE        PE 

FORMA T( 1I8,3F10.3,3F10.1,3F10.4  ) 


DZ 


TE' ,/) 


FORMAT(«  SHARP  DT  DECR E ASE ' , 2E 10. 3) 


&  £  *  £  fc  *■  fc  *;**#$:&  A:****  ###>!£:*#  **  jfc  3*  *#  #*#:£:£ 

C  'ENTIRE    CRYSTAL    CONFIGURATION' 

9690  FORMATt 14, 3F5. 2, 14) 

9691  FORMAT (9F8. 4) 
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C  'PLUCK    CONFIGURATION1 

fcA&frfcjjc&fc**  ******  »***:&:£:£:<<:#  *********** 

9690  FORMAT  (  I4,3F5.2,M  14)  ) 

9691  FORMAT( I 5,3(1X,F8.4) ,1X,I5) 

**£*#*5i;*:fc>j;*:$c#  A  ^^  A  ^*^^*  ***************  s^'t*^  *>!«**  **********  *s)t*  *  *• 


RUNTM=4*6Q-20 

START=0. 31*1 TI  ME(XX) 

DO    2    1=1,1000 

RXK( I )=0.0 

RYK{ I )=0.0 

RZK(  I  )=0.0 

VX(I )=0.0 

VY( I ) =0.0 

VZ( I)=0.0 

PKE (I ) =3.0 

PPE( I )=0.0 

PTE(  I )=0  .0 

RZKI)  =  0.0 

ISHUT=1 

NRUN=0 

READ  {  5,9010)  IH1 

READ  (  5,9023)  I H2 , DCON , AL PHA, 

READ  (  5,9030)  BU LLET , BMAS , PEXA , PEXB , I HB , THERM 

READ  (  5,9033)  TARGET , TMA S, EXA ,EXB , I HT , TEMP 

R0E2=3. 0 

READ    (     5,9350)     IHS , PL ANE , L S , I  X, I Y , I Z ,C VR , MCRO    ,DTI 
ROE=SORT (R0F2) 
ROEM    =    ROE-DTI 
R0EL2=R0EL*R0EL 
CVE=1.60E-19 
CVM=1.672E-27 
FM=1  .0F-10 
FM2=FM*FM 
CVD=CVR*1.0E-10 
CVED=CVE/CVD 
PTMAS=TMAS*CVM 
PBMAS=BMAS*CVWI 
HTMAS=0.5*PTMAS/CVE 
HBMAS=D.5*PBMAS/CVE 
VFAC=1.0 

TSAVE=BMAS/(BMAS+TMAS J 
BSAVE=TMAS/(BMAS+TMASJ 

FXA=ALOG(-EXB*CVED)+EXA 

PFXA=AL0G(-PEX8-CVED)+PEXA 

PPTC=EXP(PEXA+PEXB*ROE) 

PAC=ALCG(CVED)+PFXA 
PFPTC=EXP(PACfPEXB*ROE) 

CGD1=AL0G(DC0M) +2 . 3* AL PHA* RE 
CGD2=AL0G( 2.  0*DCON  )"+ALPHA*  RE 
CGBl=-2 .0*ALPHA*CVR 
CGB2=-ALPHA*CVR 
CGF1=AL0G(-CGB1*CVED)+CGD1 
CGF2=AL0G(-CGB2*CVED)+CGD2 

R0EA=1. 50/CVR 

ROEB=2.0/CVR 

R0EC2=R0EC*R0EC 

A(l,l)=l  .0 
A( 1,2) =ROEA 
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A( 1,3)=R0EA*R3EA 
A(l ,4)=R0EA**3 
A(  1,5)  =  EXP(EXA*-EXB*ROEA) 
A(2, 1)=1.0 
A(2f2) =ROEB 
A(  2,3)  =ROEB*ROEB 
A(2,4)=R0E3-*3 

A 1 2, 5) =EXP(CGD1+CGB1*R0EB)-EXP ( CGD2+CGB2*R0EB ) 
A(  3,  1J=0.0 
A(3, 2J=-1.0 
A(3,3) =-2. D*RCEA 
A ( 3,4) =-3. 0*ROEA*ROEA 
A(3,5)=EXP(FXA+EX8~R0EA)/CVED 
A  (  4 , 1 )  =  0  .  0 
A(4,  2)=-1.0 
A(4,3)=-2.0*R3EB 
A(4,4) =-3. D*ROEB*ROEB 

A(  4,5)=(EXP(CGF1+CGB1*R0EB)-EXP(CGF2+CGB2*R0EB) )/CVED 
CALL    CROSYM 
CP0=A(  1  ,5) 
CP1=A(2,5) 
CP2=A(3,5) 
CP3=A(4,5) 
CFO=-CPl*CVED 
CF1=-2.0*CP2*CVED 
CF2=-3.0«CP3«:VED 
5       READ    (     5,9040)     EVR , NTT ,NS ,ND , I P , I  DEE P , I  TYPE , NVAC , D1X, 
1D1Y,D1Z 
IF(NTT.EO.O)     30    TO     9999 
I0=ITYPE-1 
NTTS=NTT 
EV=EVR*1.0E+3 

*i!t)j:**^?  AAAA AAA AAAAA  ¥  <T  ***£##$  A  A*  AAA  A  ^- AAA**  **  A  A  A  A*  A*  AAA  A  A  A**** 
A** A  A* A  *AAAAA*AA*AAA*  Jr**  #^t*i{:#*A***J]( 

C  'PLUCK    CONFIGURATION' 

^•jt*^cis):if;  A  A*  A  A  AAA  AA*A**AAA*AAA*AA*A*A 

c 

C  DATA    FOR    THE    PLUCK    CRYSTAL     IS    READ    INTO    THE    COMPUTER 

C  AND    APPROPRIATE     PARAMETERS    ARE    CORRELATED.    THIS 

C  CREATES    THE    CRYSTAL    TO    BE    USED     IN    THE    SIMULATION. 
C 

READ       (5,9690)     L L , D IX, 01 Y, D 1Z , N VAC , I XNE W , I YNE W , I Z NE W 
11=60 
LL  =  I1 
IX=IXNEW 
IY=IYNEW 
IZ=IZNEW 
DO    15     1=1, LL 

RFAD       (5,9691)      I , RX ( I ) , RY ( I ) , R Z { I ) , KEEP ( I ) 
15    CONTINUE 


AAAAAAfrA^AAAAAAAAAA^cAAAAAAAAAAAAAA* 

C  'ENTIRE  CRYSTAL  CONFIGURATION' 

AA****'!c****A>f!'*A*A******AA**A*A*A*** 

C 

C  DATA  FOR  THE  ENTIRE  CRYSTAL  IS  READ  INTO  THE  COMPUTER, 
C  THIS  CREATES  THE  CRYSTAL  TO  BE  USED  IN  THE  SIMULATION, 
C 


READ   (5,9690)  LL , D1X, D1Y , D 1Z , NVAC 
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DO    15 

K=I  +  1 

J=I+2 

READ 

1RX( J), 


1=1, LL, 3 


(5,9691  )     RX(  I),  RYU  )  ,  RZ  {  I  )  ,  RX  (  K  )  ,  RY  (  K  )  ,  RZ  (  K  )  , 
RY( J),RZ( J) 


30     ILAY=IDEEP 

IFUDEEP)     35,35,40 
35    LD  =  LL 

ILAY=IY 
40    DO    45    I =1,LL 

RXK( I )=RX( I 

RYK( I J=RY( I 

RZK( I) =RZ( I 

RXI ( I)=RX( I 

RYK  I  )=RY(I 
45    RZI { I ) =RZ( I 

RXK  1)=RX(  1 

RYI (1) =RY(1 

RZI ( 1)=RZ( 1 

RXKU  )  =  RX(1 

RYKU)  =RY(1 

RZK( 1)=RZ(  1 
C 

TPOT=0. 3 

CALL    ENERGY 

BIND=-TPOT 

TE=TPOT+BIND 

TPKE=EV 

TE=TPOT+BIND+TPKE 
WRITE     (    6,9610) 
WRITF     I     6,9620)     IH2,NT 

DO    70     1=1, LL, 3 

K  =  I  +  1 

J=I+2 
70       WRITE     (6,9630)     KE EP ( I  ) , RX ( I  ) , RY ( I  )  , RZ (  I ) , PP E ( I ) , 
1KEEP(K)  ,RX(K.)  ,RY(K)  ,PZ(K)  ,PPE(K),KEEP(J)  ,RX(J),RY(J) 
1RZ(J) ,PPE( J ) 

WRITE     (     6,9643)     T PKE , TPOT ,T E, ROEL 

NPAGE=1 
WRITE     (    6,9650)    NPAGE 


C 
C 
C 
C 
C 


APPROPRIATE     IMPACT     POINTS     FOR     INTERSTITIAL     ESCAPE    ARE 
CHOSEN    AND    AN    IMPACT    POINT    GENERATOR     IS     CREATED    TO 
GENERATE    AND    THEN    TEST    POSSIBLE    DIRECTIONS    OF    ESCAPE. 


CX=3.8 

CY=0. 30 

CZ=3.6 

NTT=NTTS 

DO    3000    11=1, <t 

CX=CX+0.4 

DO    3000    JJ=1,3 

CZ=CZ+0.4 

RXSAVE=CX 

RYSAVE=CY 

RZSAVE=CZ 

AC=SORT( (CX-RXI ( 1)  )**2  +  (CY-RYI ( 1)  ) **2+ ( CZ-RZ  I  ( 1 ) )**2) 

COX=(CX-RXI (1) )/AC 

COY=CY-RYI  (1)  )  /AC 

VOL=SQRT(EV/HSMAS) 

VX(1)=V0L*C0X 

VY( 1)=V0L*C0Y 

C0Z=(CZ-RZI(1 ) )/AC 

VZU)  =VOL*COZ 

IF(NRUN.EQ.O)     GO    TO    60 
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50  DO  55  I=1»LL 
LCUT( I ) =3 
RX(  I)  =  RXI(  I ) 
RY(  I  )  =  RYI  ( I ) 
RZ( I) =RZI( I ) 
RXK(  I)=RXI  ( I  ) 
RYK{ I)=RYI (I ) 

55    RZK( I) =RZI (I) 

60    NRUN=1 

DO    65    I =2,LL 
VX(  I  )  =  0.0 
VY(  I  )  =  0.0 
VZ(I )=0.0 
PPEII J=0.0 
PKE(  I)=0.0 

65    PTE( I) =3.0 
TPOT=0.0 
NEW  =  0 

TIME=0. 0 

NT=0 

IF     (NRUN.GT.l, 

NPAGE=NPAGE+1 

95  TFAC=2.0*PTMAS 
TFACB=2.0*PBM4 
DT=1.0E-17 

100    DTOD=DT/CVD 

HDT0D=0.5*DT03 
DTOM  =  DT/PTMAS 
HDTOM=0. 5*DTCM 
DTOMB=DT/PBMAS 
HDT0MB=0.5*DTQ 

200    CALL    STEP 

IF(LCUT(1)  .GT. 
1=1 

RXKt I)=RX(  I  ) 
RYK( I }=RY( I ) 
RZK( I) =RZ( I) 
RX(  I  )=RX(I  )+DT 
RY(  I)=RY(I  )  +  DT 
RZ(I ) =RZ(I )+DT 

240    DO    245     1=2, LD 
IF {LCUT( I) .GT. 
RXK( I)=RX( I } 
RYK(  I )  =  RY(  I  ) 
RZK1I ) =RZ(I) 
RX( I ) =RX(I )+DT 
RY(  I  )  =  RY( I  )+DT 
RZ{I) =RZ(I )+DT 

245    CONTINUE 
CALL    STEP 
EMAX=3.3 
FMAX=0.0 
TIME=TIME+DT 
NT=NT+1 
IF(LCUT( 1) .GT. 
1=1 

VSS=VX( I ) 
VX( I )=VSS+HDT3 
RX(I)=RXK(I )+( 
VSS=VY( I) 
VY( I )=VSS+HDTD 
RY( I }=RYK(I )+( 
VSS  =  VZ(  I  ) 
VZ( I )=VSS+HDT3 
RZ( I) =RZK(I )+( 
PKE(I)=VX( I )*V 
EMAX=PKE(I ) 

260  FX(  I  )=0.0 
FY(  I  )  =  0.0 
FZ(I )=0.0 


0)  GO  TO  95 


*DTI*CVD 
S*DTI*CVD 


MB 

0)  GO  TO  240 


OD*(HDTOMB*FX(  I  )  +  VX( I )  ) 
OD*(HDTOMB*FY(  I  )+VY( I  )  ) 
OD-(HDTOMB*FZ(I )+VZ(  I)  ) 


3) GO  TO  245 


OD* (HDTOM*FX(I )+VX( I ) ) 
0  0-MHDTOM*FY(  I  J  +  VY(  I  )  ) 
OD*(HDTOM*FZ( I ) +V Z ( I ) ) 


0)  GO  TO  265 


MB*FX(I ) 

VX(  I  )+VSS  )*HDTOD 

MB*FY( I } 

VY(I )+VSS)*HDTOD 

MB*FZ( I ) 

VZ(I )+VSS)*HDTOD 

X(  I  }+VY(  I  )*VY( I )  +  VZ( I )*VZ(  I  ) 
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265    DO    280    1=2, LD 

IF(LCUT (  I) .GT.OJGO     TO    280 

VSS=VX( I ) 

VX( I)=VSS+HDT3M#FX( I ) 

RX< I)=RXK( I )+( VX( I )+VSS)*HDTOD 

VSS=VY( I ) 

VY(  I)=VSS+HDT3M*FY( I  ) 

RY( I)=RYK(I )+(VY( I )+VSS )*HDTOD 

VSS=VZ( I) 

VZ(  I)=VSS+HDT3M*FZ( I  ) 

RZ(I  )  =  RZK(I )  +  (VZ( I  )+VSS  )*HDTOD 

PKE( I) =VX( I )*VX( 1)+VY(I )*VY( I)+VZ( I J*VZ( I ) 

FX(  I  )  =  0.0 

FY(I )=0.0 

FZ(I )=0.0 

IF(PKE(  I  J.GT.EMAX)  EMAX  =  PKE(1) 
280  CONTINUE 

DTL=DT 

IF(EMAX.EQ.O.O )    GO    TO    285 

GO    TO    290 
285    DT=1.0E-17 

GO    TO    300 
293    DT=DTI*CVD/SQRT(EMAX1 

CTIME=0.01*ITIME(XX)-START 
300     IF (ISHUT.EO.-l )     GO    TO    400 
310     IF(NS-NT)     400,400,100 
400    CALL    PRINT 
C 

410    TPOT=0.0 

DO    450     1=1, LL 

PPE(I ) =  3.0 
450    PTE( I)=0.0 

CALL    ENERGY 

PKE(l)  =H3MAS*PKE(1  ) 

TPKE=PKE( 1) 

PTE(1 )  =  PKE(1 J+PPE(  1  ) 

DO    620    I  =2  ,LL 

PKE( I)=HTMAS*DKE( I ) 

TPKE=TPKE  +  PKE(  I  ) 
620    PTE( I) =PKE( I )+PPE( I ) 

TE=TPOT+BIND 
WRITE     (     6,9650) 

DTEST  =  (RY(  D-RYI  (1)  )**2 

IF  (DTEST.GT.  0.01)  DTEST=  0.01 
700  DO  750  1=1, LD 

DX(I )=RX(  I  J-RXI (  I  ) 

DY(  I  )  =  RY(  I  J-Rf  Id) 

DZ(I ) =RZ( I )-RZI (  I  ) 

IF    (DX( I )**2.3E.DTEST) 

IF     (DY (  I  )**2.GE.DTEST) 

IF     <DZd)**2.GE.DTEST) 

IF(PPE ( I ).GE.-3.0)     GO 

GO    TO    750 


C  'PLUCK    CONFIGURATION' 

c 

C  IN  ORDER  TO  FOLLOW  THE  LOGIC  OF  THE  DYNAMIC  PROGRAM 

C  THE  ATOMS  CREATED  BY  SUBROUTINE  PLUCK  HAVE  BEEN 

C  RENUMBERED  CONSECUTIVELY  FROM  ATOM  NUMBER  1,  BUT  THE 

C  ATOM  NUMBER  OF  EACH  ATOM  AS  IT  WAS  IN  THE  ORIGINAL 

C  CRYSTAL  HAS  BEEN  SAVED  IN  AN  ARRAY.  THE  WRITE  COMMANDS 

C  MUST  THEREFORE  PRINT  THE  ARRAY  KEEP(I)  SO  THAT 

C  PRINTED  OUTPUT  WILL  BE  IN  A  FORM  TO  ALLOW  READY 

C  COMPARISON  WITH  THE  ORIGINAL  CRYSTAL.  TO  ACCOMPLISH 


GO  TO 

720 

GO  TO 

720 

GO  TO 

720 

0  720 
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c. 
c 
c 


THIS     IN    THE    PLUCK    CONFIGURATION,    THE    STATEMENTS    NUM-      . 
BERED    720    AND    965    SHOULD    BE    CHANGED    TO    READ: 

WRITE (6,967  0)     KEEP ( I ) , DX (  I  )  , DY ( I ) , DZ{  I ) , V X (  I  ) , VY ( I  ) , 
1VZU  J  ,PKE(  I)  ,PPE(  I  )  ,PTE(I  ) 

WRITE     (6,9630)     KEEP ( I  ) , RX ( I ) , RY ( I  )  , RZ (  I ) , PPE (  I  ) , 
1KEEP(K),RX(K),RY(K),RZ(K),PPE(K),KEEP(J),RX(J),RY(J), 
1RZ(J),PPE(  J  ) 

*****************  **  *£***********:  *****  *******  **************** 


7  20 
965 


720       WRITE     (     6,9670) 
1VZ(  I),PKE(  I  ),PPE 
750    CONTINUE 

WRITE     (     6,9640) 

WRITE     (     6,9650) 

NPAGE=NPAGE+1 


I ,DX( I) ,DY(I) 
(I ) ,PTE( I) 


DZ(I) ,VX(  I)  ,VY(I  ) 


955 


965 


TPKE,TPOT,TE,ROEL 
NPAGE 


IF(NT-NTT) 
790  NS=NS+ND 

GO  TO  100 
950  CONTINUE 

NS  =  0 


790, 950,950 


WRITE 
DO  965 
K  =  I  +  1 
J=I+2 

WRITE 
1RY(K) ,R 

WRITE 

WRITE 


(  6, 
1=1, 


9620) 
LL,  3 


IH2,NT 


(  6,9630)  I, RX(  I  ),RY(  I  ) ,RZ(  I  }  ,PPE(  I  ),K,RX(K)  , 
Z(K),PPE(K) ,J,RX(J) ,RY(J),RZ(J),PPE(J) 
(  6,9640)  TPKE, TPOT,TE,ROEL 
(  6,9650)  NPAGE 


C 
C 

c 
c 
c 


AFTER  CALCULATIONS  OF  THE  DYNAMICS  ASSOCIATED  WITH 
ALL  IMPACT  POINTS  ALONG  A  (100)  DIRECTION,  THE  IMPACT 
POINT  GENERATOR  MUST  BE  INCREMENTED. 


IFCJJ.EQ.3)  CZ=3.6 
1003  IF(ISHUT)  9999,3000,3000 
3000  CONTINUE 
9999  STOP 

END 


C 
C 

c 
c 
c 
c 
c 


SUBROUTINES  B100,  PLACE,  AND  PLUCK  ARE  NOT  USED  IN 
DYNAMIC  SIMULATIONS.  SUBROUTINES  STEP,  CROYSM,  AND 
ENERGY  ARE  USED  IN  THE  SAME  FORM  AS  IN  THE  STATIC 
SIMULATIONS  AND  ARE  NOT  REPEATED  HERE.   SUBROUTINE 
PRINT  APPEARS  IN  DYNAMIC  SIMULATIONS  AS  FOLLOWS: 


C 

c 


SUBROUTINE    PRINT 


COMM ON/COM 1/RX( 1000) ,RY(1000) ,RZ(1000) ,LCUT(1000) , 
1LL,LD, ITYPE,NVAC 

C0MM0N/C0M2/IH1 ( 2 0 ) , I H2 ( 8 ) , I HS ( 10 ) , I HB ( 6 ) , I HT ( 6 ) , 
1TARGET(4) , TMA S , BULL ET( 4) , BMAS , PLANE , TEMP , THERM 

C0MM0N/C0M3/RX I (1000),RYI( 1000),RZI( 1000) ,CVR,EVR, 
1NT.TIME  ,DT,DTI  ,1  LAY 
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COMMON /COM  4/1 X, IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z 
COMMON/ COM5/ROE,ROE2,R OEM, EX A, EXB, P EXA, P E XB , FXA , P FXA , 
1I0,TSAVE,BSAVE,CX,CY,CZ 

COMMON/COM8/R0EA,ROEB,R0EC,ROEC2,CP0,CPl  ,CP2,CP3, 
lCFOtCFl, CF2.CGDl,CGD2,CGBl,CGB2tCGFl,CGF2 
C0MM0N/C0M9/RXSAVE,RYSAVE,RZSAVE 
97  10     FORMAT (40X, 10A4, /, 2  8X,  2  0A4,/) 
9720    F0RMATC9H    TAR3ET    - ,4A4 , 10HPR IMARY    -     ,4A4 ,  IX, 14HL ATT  I CE 

1    UNIT    =,F7.4,4H    ANG) 
9730    FORMAT (4X,6HMASS    = , F7 . 2  ,  13X , 6HMASS    = , F7. 2, 9X, 14HL ATTI C 
IE    TEMP     =F5.2,7H    DEG    K,,18H       THERMAL     CUTOFF    =,F5.2,3H    E 
1V/J 

9740  F0RMAT(2H     (,A^,8H)     PLANE, ,  18H       PRIMARY    ENERGY    =, 

1    F5.2,21HKEV,       CRYSTAL    SIZE     (     ,I2,3H    X     ,I2,3H    X     ,I2,3H 
1     ),,     4X,     16HVACANCY    IN     SITE     ,     14/) 

9741  F0RMAT(2H     (,A4,8H)     PLANE, ,18H       PRIMARY    ENERGY    =, 

1    F6.5,21HKEV,       CRYSTAL     SIZE     (     ,I2,3H    X     ,  12  ,3H    X     ,12, 
13H     ),,     4X,     16HINTERSTI TIAL     ( - , F 5. 2 , 2H ,- ,F 5. 2 , 2H ,+ , 
1F5.2.12H)     FROM    SITE    ,14/) 

9742  F0RMAK2H    (,A4,8H)     PLANE, ,18H       PRIMARY    ENERGY    =, 

1     F5.2,21HKEV,        CRYSTAL     SIZE     (     ,I2,3H    X    ,I2,3H    X    ,  1 2 , 3H 
1     ),,    4X,    20HR5PLACEMENT    IN    SITE     ,     14/) 

9750  FORMAT     ('     PRIMARY    START    POINT    (LU)        X=',F5.2,',       Y  =  «  , 
1F5.2,',       Z=« ,F5.2, 5X, 13, •     LAYERS    ARE    FREE     TO    MOVE', 
110X,  '10  =  ' ,12) 

9751  FORMAT     ('     OFFSET    FROM    EOUILIB     ( LU)     0X='  ,F 5. 2 , ■ , OY =' , 
1F5.2,' ,0Z=' ,F5.2,5X, 'PRIMARY    ENERGY     IS',F6.3,'    KEV',/) 

9760    FQRMAT(12H    POTENTIAL        , 6A4  ,  3X , 5HPEXA= , F9 . 5  ,2X  ,  5HP EX B=, 

1F9.5,2X,5HPFXA=, F9 . 5  ) 
9765    FCRMATt 1 2X , 6A4 , 3 X , 5HEX A    = , F9. 5 , 2X , 5HEX B    = , F9 . 5, 2X , 5HFX 

1A    =, F9.5/) 
9770     FORMATC     WHEN',F8.4,«     <    R    <',F8.4,'    THE    MATCHING    POTEN 
1TIAL    PARAMETERS    ARE',//,'  CPO    =',F10.3,',    CP1     =• 

1F10.3,',     CP2    =«,F10.3,',     CP3    =',F10.3,/,'  CFO    =' 

1E10.3,',    CF1    =»,E10.3,«,     CF2    =',E13.3,//) 
9780    FORMATC     CUT-3FF    AT',F5.2,',    WHEN    R    >    »,F6.3,'     LU ,     MOR 
1SE    POTENTIAL    PARAMETERS    ARE',     8A4,//,10X,«  CGD1    =', 

1F8.4,',     CGD2     =',F8.4,',     CGB1    =',F8.4,',     CGB2    =',F8.4, 
1«,     CGF1    =',F8.4,«,     CGF2    =«,F8-4,//) 
9790     FORMAT(10H    TIMESTEP    ,  14, 22X, 6HDT  I    =     ,     F5.4,     5H    LU, 
1,22H       ELAPSED    TIME     (SEC)     =,     E10.4,',     NEXT    TIMESTEP    IS 
1=' ,E10.4/) 

WRITE     (     6,9710)      IHS,IH1 
WRITE     (     6,9720)     TA RGET , BULL ET , CVR 
WRITE     (     6,9730)    TM AS, BMAS, TEMP , THERM 
GO    TO     (401  ,4D2  ,403,432  )  ,     ITYPE 

401  WRITE     (     6,9740)     P LANE ,E VR  ,  I  X , I Y , I Z , NVAC 
GO    TO    405 

402  WRITE     (     6,9741)     PL ANE  ,E VR , I  X , I Y , I Z , D1X , Dl Y, Dl Z, NV AC 
GO    TO    405 

403  WRITE     (    6,9742)     PL ANE , E VR ,  I  X , I Y, I Z , NVAC 

405       WRITE     (     6,9750)     RXI (  1 )  , RYI  ( 1 )  ,  RZI ( 1 ) ,  I L AY  ,  IQ 
WRITE(6,97  51)     RXSAVE,RYSAVE,RZSAVE,EVR 

WRITE     (     6,9750)      I HB, P EX  A, PEX B, P FX A 

WRITE     (     6,9765)      I HT , E XA , E X3 , F XA 

WRITE     (     6,9770)    ROEA, R0E3 , CPO , CP 1 , CP2 ,CP3 ,CF 0, CF 1 , CF2 

WRITE     (    6,9780)     R OEC , ROEB , I H2 , CGD1 , CGD2 , CGB 1 , CGB2, 
1CGF1,CGF2 

WRITE     (    6,9790)     NT  ,  DT  I  , T  IME , DT 
RETURN 
END 

BLOCK    DATA 

COMMON /COM 1/RX( 1000) ,RY( 1000) ,RZ(1000) ,LCUT(1000) , 
1LL,LD,  ITYPE, N^/AC 
DATA    R X/ 100 0*0. 0/,RY/ 1000*0. O/tRZ/1 000*0.0/, 

1LCUT/1000*0.0/ 
COMMON/ C0M3/RX I ( 1000 ) , RY I ( 1000 ) , RZ I ( 1 000 ) , C VR , E VR , 

INT, TIME ,DT,DTI  ,1  LAY 
DATA    RXI/100  0*  0.  0/,RYI/l  00  0*0.0/,  P.  ZI/100  0*0.0/ 
COMMON/ C0M6/FX ( 1 000 ) , FY ( 10  00 ) , FZ ( 1000 ) , PAC , PFPTC , FM 
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DATA    FX/ 1000*0.0/, FY/1 000* 0. 0/ t FZ/1000*0. 0 / 
END 
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